DTRA-TR-98-82-V1 


Operational Multi-Scale Environment Model with Grid 
Adaptivity, Versión 4.0,Volume 1 - Model Enhancements 


Approved for public release; distribution is unlimited. 


September 2000 



Prepared for: 

Defense Threat Reduction Agency 
45045 Aviation Drive 
Dulles, VA 20166-7517 


20010320 0/9 


DNA001-95-C-0130 
David P. Bacon, et. al. 


Prepared by: Science Applications International Corporation 

P.O. Box 1303 
McLean, VA 22102 







DESTRUCTION NOTICE: 


Destroy this report when it is no longer needed. Do not 
retum to sender. 

PLEASE NOTIFY THE DEFENSE THREAT REDUC- 
TION AGENCY, ATTN: ADM, 45045 AVIATION 
DRIVE, DULLES, VA 20166-7517, IF YOUR ADDRESS 
IS INCORRECT, IF YOU WISH IT DELETED FROM 
THE DISTRIBUTION LIST, OR IF THE ADDRESSEE IS 
NO LONGER EMPLOYED BY YOUR ORGANIZATION. 



DISTRIBUTION LIST UPDATE 


This mailer is provided to e nable DTñA to maintain current distribution lists for reporte. fWe would 
appreciate vou providina the reauested i nformation.) 

□ Add the individual listad to your distribution list 

□ Delete the cited organization/individuaL 

□ Change of address. 


Note: 

Please retum the maiting iabel from the 
document so that any additions, changas, 
corrections or deletions can be made easily. 
For distribution canceliation or more 
information cali DTRA/ADM (703) 325-1036. 


ÑAME:_ 

ORGANIZARON:_ 

OLD ADDRESS NEW ADDRESS 


TELEPHONE NUMBER: ( ) _ 

DTRA PUBLICATION NUMBER/TITLE CHANGES/DELETIONS/ADDmONS, etc.) 

(Attach Shoetttmom Space is Requimd) 


DTRA or other GOVERNMENT CONTRACT NUMBER:_:__ 

CERTIFICATION of NEED-TO-KNOW BY GOVERNMENT SPONSOR (if other than DTRA): 
SPONSORING ORGANIZATION: 

CONTRACTING OFFICER or REPRESENTATIVE: 


SIGNATURE: 













DEFENSE THREAT REDUCTION AGENCY 

ATTÑ: ADM 

45045 A VI ATI O N DRIVE 

DULLES, VA 20156-7517 


DEFENSE THREAT REDUCTION AGENCY 
ATTN: ADM 

6801 TELEGRAPH ROAD 
ALEXANDRIA, VA 22310-3398 








REPORT DOCUMENTATION PAGE 


Form Approved 
OMB No. 0704-0188 


Public reporting burden for this collection of Information is estimated to averaga 1 hour per responso, Including the time for reviewino instructions. seofcftmg exáting data sources, 
gatneting and maintainlng the data needed. and completing and revtewing the coltection of Information. Send comments regadinQ this burden estímate to any other aspect of this 
collection of Information, including suggestions for reduclng this burden, to Washington Headquarters Services, Dlrectorateffor Infonnation Operatlons and Reporta, 1215 Jefferson 
Davls Highway. Suite 1204, Arlington. VA 22202-4302. and to the Office of Management and Budget, paperworfc Reduction Project (0704-0188). Washington, DC 20503. 


1. AGENCY USE ONLY (Leave blank) I 2. REPORT DATE 


3. REPORT TYPE AND DATES COVERED 

Technicai 950615 . - 980930 


4. TULE AND SUBTiTLE 

Operational Multi-scale Environment Model with Grid 
Adaptivity, Versión 4.0, Volume 1 - Model Enhancements 


6. AUTHOR(S) 

David P. Bacon, Christopher Agritellis, Nash'at Ahmad, 
Zafer Boybeyi, Thomas J. Dunn, Mary Hall, Yi Jin, Pius Lee 
Rangarao Madala, Douglas E. Mays, Ananthakrishna Sarma, 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS( ES) 

Science Applications International Corporation 
P.O. Box 1303 
McLean, VA 22102 


5. FUNDING NUMBERS 


DNA001-95-C-0130 

62715H 

AC 

BA 

DH50464, DH52120 


8. PERFORMING ORGANIZATION 
REPORT NUMBER 


SAIC-99/1051 


9. SPONSORING /MONITOR1NG AGENCY NAME(S) AND ADDRESS( ES) 

Defense Threat Reduction Agency 
45045 Aviation Drive 
Dulles, VA 20166-7517 
CPWT / Smith 


10. SPONSORING /MONITORING 
AGENCY REPORT NUMBER 


DTRA-TR-98-82-VI 


11. SUPPLEMENTARY NOTES 

This work was sponsored by Defense Threat Reduction Agency under RDT&E RMC code B 
4662 D AC BA 50464 4400 A BA 25904D. 


12a. DISTRIBUTION/AVAILABILITY STATEMENT 


12b. DISTRIBUTION CODE 


Approved for public release; distribution is unlimited 


13. ABSTRACT (Máximum 200 words) 

This report details the status of the Operational Multi-scale Environment model 
with Grid Adaptivity (OMEGA). OMEGA is a nonhydrostatic atmospheric simulation 
system, which ineludes a graphical user interface (GUI), a grid generator, a data 
preprocessor, the atmospheric forecast model, and visualization postprocessor. It 
also ineludes terrain, land-water, and climatological surface datasets, and scripts 
that allow the model to be setup and run operationally. This report discusses the 
primary features of OMEGA with special emphasis to the improvements made under this 
contract. 

The improvements inelude a time-splitting scheme which groups the grid cells 
into masks that can use similar timesteps resulting in considerable computational 
savings; air-surface interaction routines that take into account nonhomogeneous 
surface properties; a Lagrangian partióle model that takes into account turbulent 
fluctuations in the momentum of air; and scripts which allow the automatic execution 
of the model from data ingest to results postprocessing. A major improvement to the 
system has been the dynamic adaptation capability by which the model grid adapts to 
runtime flow-characteristics to maintain the accuracy of solution and yet save 
computational óyeles. Another modification to improve runtime performance is the 
addition of Message Passing Interface (MPI) calis to implement parallelization. 


15. NUMBER OF PAGES 


16. PRICE CODE 


20. LIMITATION OF ABSTRACT 


Standard Form 298 (Rev. 2-89) 

DraOMrlKaH hu AMCI C»H 740.10 


J4.SUBJECT TERMS 

Hazard Prediction 
Unstructured Grid 
Adaptive Grid 

Atmospheric Dispersión 

Numerical Weather Prediction 

Atmospheric Prediction 

17. SECURITY CLASSIFICATION 

0F REPORT 

18. SECURITY CLASSIFICATION 

OF THIS PAGE 

19. SECURITY CLASSIFICATION 
0F ABSTRACT 

UNCLASSIFIED 

UJNCLASSIFIED 

UNCLASSIFIED 


N$N 7540-01-280-5500 


1 






















UNCLASSIFIED 

SECUfllTY CLASSIFICATION OF THIS PAGE 



UNCLASSIFIED 






Preface 


This marks the culmination of about six years of very exciting and fast-paced research. Science 
Applications International Corporation (SAIC) feels very privileged to work on the “Hazard Predication 
and Assessment Capability/Weather” project and to provide the Defense Special Weapons Agency 
(DSWA) a landmark capability to accurately predict the transport and dispersión of hazardous materials 
in air. This capability is manifested in the Operational Multiscale Environment Model with Grid 
Adaptivity (OMEGA) which significantly advances the state-of-the-art in numerical weather prediction. 
OMEGA integrates into one model the advances in the fields of computational fluid dynamics, 
numerical algorithms, atmospheric dynamics and physics, and atmospheric dispersión and transport; 
most of these advances have been the results of prior DSWA programs. The work presented in this 
report and its companion “The Operational Multiscale Environment Model with Grid Adaptivity, 
Volunte 2 - Model Validation and Verification” represents the work performed by a team forming the 
Center of Atmospheric Physics of SAIC, with the support and guidance of severa! individuáis at DSWA 
and SAIC. The SAIC’s OMEGA modeling team has enjoyed the constant support and inspiration of Dr. 
George Ullrich, Mr. Cliñon McFarland and Dr. León Wittwer of DSWA, and the expert guidance of 
Lt.Col. James Hodge, PhD., who served as the DSWA’s technical monitor for this contract and his 
successor in that role Maj. Thomas Smith. 


Table of Contents 


Section Page 

Preface iii 

Figures vii 

Tables xi 

1 Introduction 1 

2 Program Philosophy: The OMEGA Architecture 3 

2.1 Unstructured Adaptive Grids 4 

2.2 Embedded Atmospheric Dispersión Model 9 

2.3 Operational Considerations 10 

3 Solution of Goveming Equations 12 

3.1 Fundamental Equation Set 12 

3.2 Hydrodynamic Solver 14 

4 Model Physics 25 

4.1 Turbulence and the Planetary Boundary Layer 25 

4.2 Air Surface Interactions 29 

4.3 Microphysics 36 

4.4 Radiation 51 

5 Grid Generation and Adaptation 54 

5.1 Grid Definition 54 

5.2 Grid Generation 56 

5.3 Surface Characteristics Databases 61 

5.4 Static Grid Adaptation 63 

5.5 Grid Editor 66 

5.6 Dynamic Grid Adaptation 67 

6 Initial and Boundary Conditions 69 

6.1 Data Ingestión 71 

6.2 Processing of Surface Characteristics Databases 73 

6.3 Processing of Real-Time Hydrology Data 74 

6.4 PREPDAT Module 75 

6.5 GENINIT Module 77 

6.6 ANALYSIS Module 78 

6.7 Preparation of Lateral Boundary Conditions 86 


IV 




Table of Contents (Continued) 


Section 


Page 


7 Model Boundary Conditions 87 


7.1 

Lower Boundary Conditions 

87 

7.2 

Lateral Boundary Conditions 

88 

7.3 

Upper Boundary Conditions 

89 

Atmospheric Dispersión Model Structure 

91 

8.1 

OverView 

91 

8.2 

Eulerian Model 

91 

8.3 

Lagrangian Partióle Model 

91 

8.4 

Prohabilistic Puff Model 

96 

8.5 

Lagrangian Model Modes of Operation 

100 

Dynamic Grid Adaptation 

101 

9.1 

Introduction 

101 

9.2 

Methodology 

102 

9.3 

Adaptation Criteria 

102 

9.4 

Pseudo-Laplacian Interpolation 

103 

9.5 

Adaptation to Atmospheric Dispersión 

104 

9.6 

Effects of Improving the Underlying Terrain 

105 

9.7 

Adaptation to Velocity Deformation 

108 

Parallelization 

111 

XOMEGA Graphical User ínterface 

117 

11.1 

OverView 

117 

11.2 

OMEGA Shell Script Interface 

118 

11.3 

XOMEGA 

120 

11.4 

Set Time Domain 

126 

11.5 

Set Spatial Domain 

128 

11.6 

Set ADM Sources 

134 

11.7 

Run Grid Generator 

138 

11.8 

Set Meteorology 

140 

11.9 

Fetch Meteorology 

143 

11.10 

Run Preprocessor 

145 

11.11 

Run OMEGA Model 

147 


v 



Table of Contents (Continued) 


Section Page 

12 OMEGA Output and Visualization 150 

12.1 OMEGA Output 150 

12.2 X-Grid - OMEGA Visualization Tool 150 

12.3 OMEGA Simulation Image Maker (OSM) 154 

13 Conclusions 160 

14 References 161 

Distribution List DL-1 


vi 


Figures 


Figure Paa e 

o 

1 OMEGA grid element 5 

2 OMEGA coordínate system and vertical alignment of OMEGA grid 5 

3 Vertex addition and reconnection 7 

4 Vertex deletion and reconnection 7 

5 Vertex relaxation and edge bifurcation 8 

6 OMEGA represents a balance of physical, numerical, and operational 

considerations 11 

7 Arrangement of the variables on a 2-dimensional unstructured-grid. The generic 

advected variable, vj/, is placed at cell centroids. The subscripts L and R 
designates their left-hand-side and right-hand-side cell placements. The velocity 
is decomposed into normal and parallel components (not to scale). The flux F is 
co-located with V n at the face centers, and points only in the normal direction 14 

8 Rotating cone advection test 17 

9 Rotating cone advection test after five rotations 18 

10 Error scaling with grid resolution. Three measurements of deviation from the 
original cone are used: root-mean-square error over the whole mesh, root-mean- 
square error over an area extended by the original cone, and deviation in the 

height of the cone 19 

11 The ghost triangle used in the calculation of the Laplacian 21 

11 A local coordínate system defines the horizontal plañe passing through the 

centroid of each triangle. 21 

13 Advection test consisting of a passive scalar in a steady-state flow field 
composed of six counter-rotating vórtices (left). The analytic solution at a 
given scaled time (from Staniforth) is shown (center); along with the second 

order solution at the for the OMEGA advection solver (right) 23 

14 Coarse and two finer sub-domains for testing time splitting technique 24 

15 Initial and final positions of the cone 24 

vii 



Figures (Continued) 


Page 

16 Microphysical processes parameterized in OMEGA 37 

17 Using an appropriately chosen set of criteriz, OMEGA can generate a grid 

that captures the important physical features 55 

18- Initial rectangular grid 57 

19 Grid refinement 58 

20 Grid after reconnection. The deleted edges are shown as dashed lines 

and the new edges are shown as thin lines 59 


21 The vertical levels of the OMEGA grid are constructed by extending 
radiáis through each of the surface vértices. A set of stretchable layers 


exists near the surface with constant thickness layers above 60 

22 Dataseis used to initialized OMEGA 64 

23 OMEGA grid for the Middle-East región used for gulf war analysis. The 

grid has been refined around data sources 65 

24 Refined to regional interim emissions inventory 66 

25 Test of the dynamic adaptation routine. The adaptation was to the presence 
of test particles (shown in red) initialized in the Southwest comer of the 

grid and traveling on a prescribed trajectory 68 

26 Data flow between the preprocessor and other OMEGA components 70 


27 A particle Al, located in one element of the OMEGA unstructured grid, 
is evaluated using A x B. The cross product between the vector formed 
from the first vertex location to the current particle location, A , and the 
first vertex location to the second vertex location, B indicates a positive 

valué of area with respect to B ■ 94 

28 A particle, XI, located in one element of the OMEGA unstructured grid, 
is evaluated using A x B. The cross product between the vector formed from 
the first vertex location to the current location to the current particle location, 

A , and the vector formed from the first vertex location to the second vertex 
location, B indicates a negative valué of area with respect xo B . 


vi i i 


94 



Figures (Continued) 


Figure Page 

29 Computational grid for Case A simulation (left) and Case B simulation (right) 105 

30 Partióle locations for Case A simulation (left) and Case B simulation (right) 106 

31 Time evolution of pressure perturbaron field for the growing mountain case. 

Vertical cross sections at initialization (left), 8 minutes (center), and 15 

minutes (right) 107 

32 Pressure perturbation fields for lOOOm static mountain case (left) and 350m- 

lOOOm growing mountain case (right) after 30 minutes into the simulation 107 

33 Particle locations for 350m static mountain simulation (left); lOOOm static 

simulation (center) and 350-lOOOm growing mountain simulation (right) 108 

34 Simulations for the Four Comers, USA. Static (Case A) simulation results 
after 12 hours into the simulation are shown in the left panel and dynamically 
adapting simulation (Case B) results after 12 hours into the simulation are 

shown in the right panel 109 

35 Adaptation to velocity deformation. Initial grid (left) and sea level pressure 

(right) at the start of the simulation 110 

36 Adaptation to velocity deformation. Final grid (left) and sea level pressure 

(right) after 12 hours into the simulation 110 

37 Speedup projection of the OMEGA code by a power Fortran compiler 112 

38 File structure for the parallelized OMEGA code 113 

39 The kl Om model with one domain, 9 sub-domains and 16 sub-domains 114 

40 Run time plotted as a function of processor allocation 115 

41 Run time and parallel efficiency plotted as a function of processor allocation 115 

42 Run time and parallel efficiency plotted as a function of processor allocation 

for a mainframe and a workstation cluster 116 

43 SelectCase 119 

44 ViewCase 119 


IX 



Figures (Continued) 


Figure 


Page 

45 

XOMEGA 

121 

46 

Set preferences dialog box 

122 

47 

Map operations dialog box 

124 

48 

Draw area features 

125 

49 

Set time domain dialog box 

127 

50 

Set spatial domain dialog box 

129 

51 

Set refinement disks dialog box 

133 

52 

Set ADM sources dialog box 

135 

53 

Run grid generator dialog box 

139 

54 

Set meteorology dialog box 

142 

55 

Fetch meteorology dialog box 

144 

56 

Run preprocessor dialog box 

146 

57 

Run model dialog box 

148 

58 

XGRID terrain view upon initial data load 

151 

59 

Display of map and contour options 

152 

60 

Zoomed in, with wind arrows 

152 

61 

XGRID vertical cross section 

153 

62 

XGRID profile 

153 

63 

XGRID skewT-logP diagram 

155 

64 

Field of regard and particle trajectories 

156 

65 

Concentration valúes computed from the ADM Lagrangian tracers 

156 


X 




Figures (Continued) 


Figure 


Page 

66 

OSIM stage one 

157 

67 

OSIM stage two 

158 

68 

OSIM stage three 

159 


XI 



Tables 


Table Page 

1 Soil parameters of the USDA (United States Department of Agriculture) texturai classes 34 

2 High-resolution OMEGA terrain and land/water databases 56 

3 Dataseis used to determine surface characteristics 61 

4 Summary of NCEP GRIB and BUFR data files ingested by PREPDAT 73 

5 Summary ofFNMOC data files ingested by PREPDAT 73 

6 Summary of data types which can be ingested by PREPDAT 76 

7 Observation increment error thresholds 80 

8 Recommended valúes for user specified OI variables 83 


xii 



Section 1 


Introduction 


The Operational Multiscale Environment model with Grid Adaptivity (OMEGA) and its embedded 
Atmospheric Dispersión Model (ADM) was designed to support critical studies and forecasts of 
hazardous atmospheric dispersal. OMEGA is a multi-scale, non-hydrostatic atmospheric simulation 
model with an adaptive grid that permits a variable spatial resolution that ranges from 100 km to 1 km 
without wave-reflecting intemal bonndaries. The model contains a 2.5 level explicit boundary layer 
formulation, surface layer physics with múltiple soil layers and 12 soil types, a sophisticated explicit 
microphysics package with five water species and a four dimensional data analysis scheme based on the 
optimum interpolation technique. ADM is an atmospheric dispersión model that computes the 
dispersión of aerosols and gases. This report is split into two volumes. Volume 1 contains the 
OMEGA/ADM atmospheric forecast and dispersión system. Model validations against analytic and 
idealized test problems, and against meteorological and dispersión dataseis are described in Volume 2. 

This report is intended to replace the documentaron of OMEGA v. 1.0 (Bacon, et al, 1996). This report 
documents OMEGA v. 4.0, which ineludes capabilities for dynamic grid adaptation and a versión, which 
has been parallelized via domain decomposition using the Message Passing Interface (MPI). Since the 
first report, many modifications, which are listed below, have been made to the model and its 
operational system. The major changes in the OMEGA system are: 

• Múltiple layer cell masking to improve model execution speed, 

• Automatic incorporation of heterogeneous surface properties from archived databases, 

• Improved boundary conditions (top, bottom, and lateral), 

• Expanded data pre-processor options to permit a wider variety of input data, 

• Splitting and merging of puffs and the development of a probabilistic transport model, 

• Improved and expanded user interface, 

• Development of Scripts to support hands-free automated daily model execution including 
automatic generation of a standard graphics package, 

• Explicit boundary layer formulation using a 2.5 level K-E scheme, 

• Improved surface energy physics for improving surface forecasts with two soil layers 
and a canopy layer. Twelve soil types are allowed in the model, 

• Inclusión of microphysics allowing five species of water and a cumulus 
parameterization for coarser resolution areas, 

• Dynamic grid adaptation to increase resolution wherever and whenever needed, 

• A faster time-split solver to reduce the OMEGA run time, 

• User fiiendly data ingest routines, 

• XGRID system visualization tool, 

• Packed Binary data format to reduce the storage requirement without sacrificing data 
accuracy, and, 

• Model parallelization to run on multi-processor Computer systems. 
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The basic OMEGA grid architecture is described in Section 2. Section 3 describes goveming equations 
and hydrodynamics solvers. Model physics is presented in Section 4. The OMEGA grid generator and 
dynamic adaptation for heterogeneous surface properties are presented in Section 5. Sections 6 and 7 
describe model boundary conditions, data quality control methods and optimum interpolation data 
analysis technique. The atmospheric dispersión model used for particle and plume dispersión is 
described in Section 8. Dynamic adaptation for moving meteorological systems or for dispersing 
plumes is described in Section 9. Section 10 describes our recent progress in parallelizing the OMEGA 
code for massively parallel computers. A detailed description of the user friendly XOMEGA Graphical 
User Interface (GUI) is presented in Section 11. Section 12 describes the XGRID visualization package 
for OMEGA/ADM. It also contains a machine independent data storage technique (Packed Binary 
Format) developed for OMEGA; this technique substantially reduces data storage requirements for the 
system. As mentioned above, Volume 2 of this report contains model validation against analytic and 
idealized test problems, and against meteorological and dispersión data sets. 
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Section 2 


Program Philosophy: The OMEGA Architecture 


In order to improve the fidelity of hazardous transpon models, it is essential that the meteorological 
forecast itself be improved. This is because the modeling of atmospheric dispersión involves virtually 
all scales of atmospheric motion from microscale turbulence to planetary scale waves. The current 
operational atmospheric simulation systems (Hoke, et al., 1989; Janjic, 1990; Mesinger, et al., 1988) are 
scale specific and cannot resolve the full spectmm required for the accurate forecast of local scale 
phenomena. Even with recent advances in computational power (McPherson, 1991), the current 
architecture and physics of today’s generation of atmospheric models cannot fully simúlate the scale 
interaction of the atmosphere. Recently, several groups have started the development of non-hydrostatic, 
nested (multiply nested in some cases) atmospheric models (Dudia, 1993; Skamarock and Klemp, 1992); 
however these represent an incremental evolutionary path in atmospheric simulation. 

A capability such as OMEGA/ADM has been desirable for many years. Had it existed, it would have 
found broad use in the Defense Nuclear Agency (DNA) Global Effects Program that was mounted to 
study Nuclear Winter related issues. It would also have proven useful in forecasting the trajectories of 
the dust and ash injected into the atmosphere by numerous recent volcanoes including Mt. St. Helens, El 
Chichón, Mt. Pinatubo, Mt. Redoubt, and Mt. Spur. OMEGA/ADM did not exist primarily because of a 
perception that there is no present and continuing need for a permanent forecast capability. This 
. -perception changed during OPERATION DESERT SHIELD/STORM (hereafter referred to as DESERT 
STORM). Science Applications International Corporation (SAIC) proposed the creation of 
OMEGA/ADM as a method for providing hazardous aerosol and gas dispersión prediction support for 
future US military operations. 

OMEGA/ADM was conceived as an operational tool to support US forces in the field. It was to form 
the kemel of an Operations Center. As in DESERT STORM, the Operations Center would be centrally 
located and would receive worldwide meteorological information. As will be seen, however, the rapid 
increase in computational power available in high performance workstations has led to the concept of 
using OMEGA in the field as well, albeit at a degraded performance level from that available from 
supercomputers. 

The major advantages of OMEGA/ADM over the current state-of-the-art inelude the ability to resolve 
the surface terrain down to scales of 1 km and along with that the local perturbations on the larger scale 
wind field. This local wind field perturbation is of extreme importance in determining the trajectory of 
an aerosol release. In order to calcúlate this local perturbation, however, it is important to inelude all of 
the physical parameters and processes that affect the local flow. These inelude not only the topography, 
but the land use, the land/water composition, the vegetation, the soil moisture, the snow cover (if 
appropriate), and the surface moisture and energy budgets. The inclusión of this additional physics, 
some of, which is only appropriate because of the increased spatial resolution, represents an additional 
advance in the state-of-the-art. 
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OMEGA, with its embedded ADM was conceived to advance the state-of-the-art in predicting the 
transport and diffusion of hazardous releases. The bulk of hazardous releases occur near the surface, are 
dispersed primarily in the PBL, and are strongly influenced by surface features (Sherman, 1978; Paegle, 
et al., 1984). These hazardous releases often require emergency response. Effective emergency 
response, in tum, requires the highest possible resolution of both the atmospheric State as well as the 
aerosol concentration. OMEGA is based upon an adaptive unstructured grid technique (AGARD, 1992) 
that makes possible a continuously varying horizontal grid resolution ranging from 100 km down to 1 
km and a vertical resolution from a few tens of meters in the boundary layer to 1 km in the free 
atmosphere. This feature allows one to obtain the highest possible resolution of the atmosphere as well 
as the hazardous aerosol/gas concentration. 

The basic philosophy of the OMEGA/ADM model development has been the creation of an operational 
tool for real-time hazard prediction. The model development has been guided by two basic design 
considerations in order to meet the operational requirements: (1) the application of an unstructured mesh 
numerical technique to atmospheric simulation; and (2) the use of an embedded atmospheric dispersión 
algorithm. For many years, SAIC has been developing new algorithms and techniques in computational 
fluid dynamics (CFD); SAIC has also been studying cloud scale and mesoscale dynamics, 
thermodynamics, and microphysics, and atmospheric transport at cloud scale and mesoscale. 

OMEGA/ADM represents the marriage of these efforts. 

The use of an unstructured mesh in OMEGA permits the use of variable grid resolution in a natural 
fashion, as well as providing a relatively straightforward basis for the introduction of dynamic grid 
adaptivity. The use of an embedded dispersión model eliminates the input/output (I/O) restrictions that 
would otherwise domínate the operational System at 1 km resolution. In this section, we discuss first the 
unique unstructured grid used in OMEGA. We then discuss the important considerations that led to the 
development of ADM. Finally, we discuss the necessary balance between physical fidelity, numerical 
accuracy, and operational constraints that affected the overall design of the OMEGA system. 

2.1 Unstructured Adaptive Grids. 

We will discuss the OMEGA model grid structure in more detail in Section 3; however, it is so 
intricately enmeshed in the overall design that it is important to set the stage for our discussion here. 

The flexibility of unstructured grids facilitates the gridding of arbitrary surfaces and volumes in three 
dimensions. In particular, unstructured grid cells in the horizontal dimensión can increase local 
resolution to better capture topography or the important physical features of atmospheric circulation and 
cloud dynamics. The underlying mathematics and numerical implementation of unstructured adaptive 
grid techniques have been evolving rapidly, and in many fields of application there is recognition that 
these methods are more efficient and accurate than the structured logical grid approach used in more 
traditional codes (Baum and Lohner, 1993; Schnack, et al., 1993). To date, however unstructured grids 
and grid adaptivity have not been used in the atmospheric Science community (Skamarock and Klemp, 
1992). Atmospheric model problems characterized by large model domains, long time integration, and 
operational time constraints have not been viewed as viable candidates for the emerging CFD grid 
technologies. OMEGA represents the first attempt to use this advanced CFD technique for atmospheric 
simulation. 
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Figure 2. OMEGA coordínate system and vertical 
alignment of OMEGA grid. 
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OMEGA is based on an unstructured, triangular prism (Figure 1) computational mesh. This mesh is 
unstructured in the horizontal dimensión and structured in the vertical dimensión (Figure 2). The 
rationale for this mesh is the physical reality that the atmosphere is highly variable horizontally, but 
always stratified vertically. While completely unstructured three-dimensional meshes have been used 
for other purposes (Baum, et al., 1993; Luo, et al., 1994), the benefít of having a structured vertical 
dimensión is a significant reduction in the computational requirements of the model. 

In duscussing unstructured grids, it is necessary to define the elemental objects that describe the 
properties of a volumetric mesh. The lowest order object of a grid is the vertex, which is specified by its 
position (x,y,z). An edge conveys the connectivity of the mesh and is defined by the indices of its 
starting and terminating vértices. The face represents the interface area between adj acent volumetric 
cells and can be described from the list of edges that bound it, or by the sequence of vértices that form its 
extremities. The cell or control volume is in tum specified by the list of faces that contain it. 

On an unstructured grid, the number of edges that meet at a vertex is arbitrary. Consequently there is no 
longer a simple algebraic construct that can be used to deduce the relationship of indices for the various 
elemental objects, as in the case of structured grids that have been used to date as the basic structure for 
atmospheric and ocean circulation models. The formation of the grid is tied to the actual solution of the 
model equations and to the topography. This means that the initial grid can be readily adapted to surface 
features or other fixed terrain features as well as the initial weather. 

An important adjunct to the unstructured triangular grid methodology is the requirement that the model 
has to compute the normal to each face in order to calcúlate the fluxes across the boundaries. Since 
these normáis must be computed, there is no benefit from orienting the grid in any particular fashion, so 
long as the numerical resolution of the hardware is sufficient to evalúate the critical fluxes. (We shall 
discuss this in more detail later in this section.) This leads to a natural separation between the coordínate 
system for the fundamental equation set and the grid structure. The coordínate system can be as simple 
as possible (such as Cartesian) while the grid structure, in this coordínate system, is extremely complex. 
OMEGA uses a rotating Cartesian coordínate system, but the grid structure is terrain following. Figure 2 
shows a rotating Cartesian coordínate system in which the origin is the center of the Earth, the z-axis 
passes through the North Pole, the x-axis passes through the intersection of the Equator and the Prime 
Meridian, and the y-axis is orthogonal to both. 

In this coordínate frame, the equations of motion are in their simplest possible form (without going into 
a non-rotating frame that would lead to unusual boundary conditions as the surface terrain moved 
through the grid) with only two temas that are somewhat non-conventional: gravity and the Coriolis 
acceleration. Gravity in this frame is directed in the — r direction (r representing the unit vector along 
an outward radial that implies that it potentially has components in all three coordinate directions. The 
Coriolis forcé is by definition 2 pQ x v and likewise has components in all three coordinate directions. 

An important aspect of the OMEGA grid structure is the fact that all the cells along a radial possess the 
same footprint. This is accomplished by creating a surface grid and then projecting radiáis from the 
center of the Earth through all of the vértices. Each layer is then created by specifying a set of vértices 
along each radial (Figure 2) (the number of vértices along any radial are the same). 
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Figure 3. Vertex addition and reconnection. 

Because the OMEGA grid structure results in the mixing of the Earth-relative horizontal and vertical 
components, it is essential that the numerical scheme be able to sepárate these. Given a grid structure 
that may be a few meters in vertical resolution and a few kilometers or a few tens of kilometers in 
horizontal resolution, the numerical resolution must be accurate to better than one (1) part in 10 5 . 

The adaptation of an unstructured grid takes place through a variety of grid operations. The first is 
vertex addition and is usually followed by a vertex reconnection step. Figure 3 illustrates these two steps 
when some activity that would indicate a need for more resolution is noted in two cells. The vertex 
addition step is accomplished by adding a vertex at the centroid of each affected cell and connecting it to 
the vértices of the cell. The reconnection step then involves the evaluation of each new cell to see if it is 
possible to create grid cells with lower aspect ratio by removing an edge and reconnecting the altemative 
vértices. 

Figure 4 shows the reverse process in which the grid is coarsened through the process of vertex deletion. 
This is also usually followed by a vertex reconnection step. It is important to note that even though the 
grid adaptation routines may create an apparent motion of the grid, it does not, in fact, move; rather the 
goal is to refine the grid in advance of any important physical process that could require additional grid 



Figure 4. 


Vertex deletion and reconnection. 









resolution, and to coarsen the grid behind the región. This differentiates the method from the adaptation 
techniques described by Dietachmeyer and Drogemeier (1991) that used vertex movement to adapt to 
atmospheric features. 

A different type of process is shown in Figure 5. In this figure we show vertex relaxation, in which the 
vértices are allowed to move as a mass-spring system, and edge bifurcation that is equivalent to vertex 
addition in the special case of a boundary cell. 

Currently, the OMEGA grid automatically adapts to static features of the underlying topography and 
land/water boundary. In addition, other factors (either related to the surface or the atmosphere) can be 
included through a user-defined field. The code can also be set to adapt the grid during runtime. This 
enables atmospheric features that require additional grid points for adequate simulation to be resolved by 
a grid that refines or coarsens based on a given set of adaptation rules. In Section 3 we will describe the 
procedures for refinement (vertex addition), coarsening (vertex deletion), and reconnection of OMEGA 
grid system. Each of these procedures is applied to generate the best possible grid based on 
considerations of resolution and element shape. 

The adaptation procedures are the essential features that give unstructured grid algorithms for partial 
differential equations their great power. This happens on two different levels. First, the total number of 
grid points necessary to perform a successful numerical computation that recovers the correct physics 
can be greatly reduced. By this we mean that the recovery of the model physics at the smallest scale 
length resolved does not require the complete domain to have the same resolution. The resources of the 
numerical and computational machinery are focused on the regions of importance. This saving is 
especially true for three dimensional hydrodynamic problems where the resulting economy, in our 
experience has made the difference between tractability and intractability. The grid can be coarse where 
the circulation is regular and smooth, but greatly enhanced where there are sharp gradients, where there 
are discontinuities in the flow, where topographic features are important, or where model physics or 
dispersión source terms require fine resolution. 




Figure 5. Vertex relaxation and edge bifurcation. 
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Second, the flexibility of unstructured grids to adapt to the transient physical solution improves the 
fidelity of the finite difference or finite element numerical schemes, eliminating some of the errors that 
plague existing models. The improvement comes from aligning the grid structure to the flow and from 
the local refmement of the grid in the vicinity of rapidly changing horizontal spatial structures in the 
atmosphere. 

2.2 Embedded Atmospheric Dispersión Model. 

Over the past decade, there has been an increasing interest in using detailed meteorological information 
in aerosol and gas transport and diffusion models to improve our understanding of atmospheric 
dispersión. The detailed meteorological information in such models could be obtained either from 
meteorological observations or from meteorological models. Until recently, the conventional weather 
observing network has been used to obtain the necessary meteorology for dispersión models (c.f, Pack 
et al., 1978; Artz, et al, 1985). These and other studies, however, have demonstrated that the weather 
observing network is too coarse in both time and space to accurately depict local circulations. 

An altemative approach to overeóme the limitations of the existing weather observation network is to 
use atmospheric conditions predicted by meteorological models. Therefore, advanced atmospheric 
dispersión models have involved the coupling of the dispersión model with a meteorological model. 

One such coupling was accomplished by McNider et al. (1988), who studied the influence of diumal 
variation and inertial boundary layer oscillations on long-range dispersión. Pitts and Lyons (1992) 
applied their coupled model to the urban area of Perth, Australia. Still another coupling was performed 
by Boybeyi et al. (1995) with an application to the Bhopal gas accident in India. 

Although this coupling approach has been sufficient for these and many other studies up until now, the 
major problem with this approach is that atmospheric dispersión models can only use the results from 
the meteorological models at some spatial and temporal resolution. As the grid resolution is increased, 
the data requirements of this approach grow as the inverse cube of the grid resolution. This can be seen 
by recognizing that the time between “snapshots” of the wind field must be less than the transit time 
across a cell. Thus, if the spatial resolution increased by a factor of two, then the data must be provided 
twice as frequently. An increase in the grid resolution of a factor of two, therefore results in an increase 
in the data requirements for the model of a factor of eight. Many current aerosol transport and diffusion 
models use analyzed gridded wind fields with roughly 100 km grid resolution that are output every 
12 hours. To improve the aerosol transport capability by increasing the grid resolution to 1 km would 
increase the data handling requirements of these models by a factor of 10 6 . 

Models-3, the next generation air quality model of the US Environmental Protection Agency (EPA) 
(Dennis and Novak, 1991) is an example of the current coupling paradigm carried out to higher 
resolution. The typical performance characteristics of an EPA study give this approach some merit. The 
meteorological component of the modeling system is re-used many times over in these studies. A typical 
study may vaiy pollutant emissions assumptions and alter the atmospheric chemistry package for fixed 
atmospheric conditions in order to develop an understanding of the driving factors in that particular air 
quality situation and to aid in the development of control strategies. Under this paradigm, however, first 
a meteorological model is used to simúlate the atmospheric conditions that are output at selected times, 
and then a chemistry and transport model evaluates the fate of the pollutants under these atmospheric 
conditions. The Models-3 program uses the NCAR/Penn State University MM5 model 
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coupled to the Atmospheric Chemistry Model (a suite of chemistry and transport models); other similar 
Systems inelude the Colorado State University Regional Atmospheric Modeling System (RAMS) 
coupled to a particle dispersión model (McNider, et al., 1988; Pielke et al., 1991). 

A major drawback to this approach, however, is the large amount of data that must be output by the 
meteorological model and ingested by the particle dispersión model. As higher and higher resolution is 
implemented, this method quickly becomes bound by its input/output requirements. OMEGA 
overcomes this problem by embedding its atmospheric dispersión model (ADM) into the meteorological 
model. OMEGA/ADM has both Eulerian (grid based) and Lagrangian (grid free) particle and puff 
transport and diffusion algorithms. These algorithms have the benefit of the full OMEGA simulation 
fields (temperature, moisture, wind, turbulence, etc.) at each timestep. In addition, the flexible grid 
structure of OMEGA permits the use of high resolution where it is naturally needed (e.g., in regions of 
complex terrain, near land/water boundaries, near emission sources) without burdening the model with 
high resolution everywhere. We will discuss the ADM in better detail in Section 8. 

2.3 Operational Considerations. 

The ultímate consideration in the OMEGA design has always been world-wide operational utility. 
OMEGA has been designed with a delicate balance between physical fidelity, numerical accuracy, and 
operational constraints (Figure 6). For example, the design point of 1 km for the highest grid resolution 
was based upon a recognition that although world-wide terrain elevation and land/water boundary 
information is available at even higher resolution, world-wide albedo, land use, vegetation, soil texture, 
and soil moisture, as well as the world wide meteorological data observations, are rarely available at 
even 1 km resolution. 

The world-wide real-time operational requirement influenced a number of design decisions. Real-time 
operational use for a reasonable forecast period dictated a realistic máximum resolution of 1 km. This 
determined one level of requirements for the physics that had to be contained in OMEGA. (We discuss 
the OMEGA physics in more detail in Section 4.) World-wide utility required that any physical models 
that were incorporated into OMEGA would have to have world-wide applicability and could not rely on 
datasets that were not available with this coverage. This, for example, eliminated the use of a second- 
order closure turbulence model since the Reynolds stress is not operationally available for either initial 
or boundary conditions. 

As part of the real-time operational requirement, SAIC recognized that a great deal of automation was 
desirable in the OMEGA system. This resulted in the creation of a highly automated grid generator 
(discussed in Section 5), an automated meteorological and surface data assimilation system (discussed in 
Sections 6 and 7), and a user friendly X-windows and Motif™ based (GUI) (Section 11) and graphic 
post processors (Section 12). 

The world-wide requirement forced SAIC to develop a number of world-wide databases including 
roughly 1 km resolution terrain elevation, standard deviation of terrain elevation (an indicator of surface 
roughness), and water fraction, coarser 5 arc-minute databases of the same valúes, and datasets for land 
use, vegetation, and soil texture at various resolutions. These are discussed in Section 5. 
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Figure 6. OMEGA represents a balance of physical, numberical and operational 
considerations. 

With this as a basis, the design objective of OMEGA is a tool for real-time hazardous prediction capable 
of performing a full coupled atmospheric forecast and aerosol or gas dispersión calculation at 10 times 
real-time. The acceptance level is half of this, or 5 times real-time. As will be seen, the tremendous 
flexibility of the OMEGA modeling system permits the model to meet these constraints on platforms 
ranging from CRAY supercomputers to IBM R/S 6000 workstations - albeit by trading of resolution or 
size of the computational domain. The remainder of this document provides the details of the OMEGA 
model, its physics, numerics, and the operational design of the system. 
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Section 3 

Solution of Governing Equations 


3.1 Fundamental Equation Set. 

In this section, we document the fully elastic non-hydrostatic equation set used in OMEGA, including 
the applicable assumptions. For brevity, we classify the five mixing ratios for water substances into two 
groups: precipitating water substances, Q p (where the subscript p designates either rain, snow, or hail), 
and nonprecipitating water substances, Q„ (where n designates either ice crystals, water vapor, or water 
droplets). Furthermore, we cast the equations in their conservative form consistent with the fully elastic 

mass-conservation equation (dpj dt = —V • p \). This form is better suited for the state-of-the-art 
upwind advection schemes we have used, which have significantly less numerical diffusion (important 
for the finer resolution we are trying to achieve) than the schemes used in other non-hydrostatic models 
such as the Terminal Area Simulation System (TASS) (Proctor, 19S7). 

We begin by decomposing the atmospheric pressure and density into a time-invariant hydrostatic base 
State and a perturbation upon that State such that 

p(x, y, z, t) = p 0 (x, y, z ) + p (x, y, z,t ) (1 a) 

and p(x,y,z,t) = p 0 (x,y,z)+p'(x,y,z,t) (Ib) 

with 

x,y,z) = -gp 0 (x,y,z ). (2) 

óz 

Assuming an Eulerian grid, we then have the following equation set for the conservation of mass, 
momentum, and energy. 


Conservation of Mass: 

& + ?.(pX) + F e = 0 

^ L +9-(pe.v)=pM.+F a 

d ( 3 ) 

- e § l +'(pQp)= pM r +|(e,w„p)+ F Q¡ 

^+v(pe„v)=pM„ +|■(e.M'.pJ+Fa 

In equations 3.3, p represents the dry air density, calculated from the total density as: 

P = P,o'a,/{l- + Q V a P cr)- ( 4 ) 
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Conservation of Momentum: 


^ + V-(pVv)=-Vp'-(p-p 0 ) g r + 2pñxV + F„ (5) 

where the density used in the calculation of momentum is the total density. 

Conservation of Energy: 

^+V-(£V)=^S(¿A) +S « +F . (6) 

Terms have been arranged such that the conservative advection terms appear on the left side of each 
equation. The source terms on the right side of the momentum equation also inelude buoyaney and 
gravitational effeets, -(p-p 0 )gr (where r is the radial unit vector), and the Coriolis forcé, (2 pQ. x Y). 

Sub-grid scale turbulence contributions, F, are discussed in detail in Section 4.1. For the remaining 
equations, T is the temperature, L¡ and S¡ denote the latent heat and rate of phase conversión of either 
vaporization, fusión, or sublimation, and W p represents the terminal velocity of each of the precipitating 
water substances. S¡ depends on the microphysics that govems the rate of phase transitions and W a 
depends on the assumed size distribution and mass of the hydrometeors. M n and M p are the non- 
precipitating and precipitating microphysics source terms (see Section 4.2). Finally, subscript a refers to 
transport of aerosol or gas. 

The initial energy density E is calculated from the total pressure as: 



and throughout the simulation pressure is diagnosed from energy density using this relationship. 
Potential temperature, 0, and temperature, T, are related by Poisson’s equation. 


where 



( 8 ) 


p ref = lOOOmb, 

K = — = 0.286, 

c„ 


(9) 


R d being the gas constant for dry air and c p is the specific heat of dry air at constant pressure. Finally, 

Sr represents the contribution of radiation flux to heating the atmosphere and will be described in 
Section 4.4. 
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3.2 Hydrodynamic Solver. 

The hydrodynamic elements of the OMEGA model are based on numerical methods of solution of the 
Navier-Stokes equations on an unstructured grid in the horizontal direction and a structured grid in the 
vertical. We use a standard split-operator methodology, calculating advection terms explicitly and 
diffusion terms implicitly. In the calculation of momentum, the pressure gradient, Coriolis, and 
buoyancy terms are calculated explicitly along with the advection terms. An implicit vertical filter and 
an explicit horizontal filter are applied to the vertical momentum. The calculation of the new 
momentum at each time-step thus involves several steps, which are described below. All implicit 
operations are performed by tri-diagonal matrix inversión. 

3.2.1 Advection - Implementation of the Finite Volume Upwind Solver. 

In this section we describe the implementation on unstructured grids of the multidimensional positive 
definite advection transpon algorithm (MPDATA) originally developed on regular grids by 
Smolarkiewicz (1984), Smolarkiewicz and Clark (1986), and Smolarkiewicz and Grabowski (1990). 

The resulting scheme is second-order-accurate in time and space, conservative, combines the virtues of 
the MPDATA (e.g., ability to separately ensure monotonicity and positive definiteness) with the 
flexibility of unstructured grids (Baum and Lohner, 1989, Peraire et al., 1990), is compatible with 
adaptive mesh algorithms (Fritts, 1988), and can run efficiently on highly parallel computers. Below, 
we describe the essential methodology and demónstrate the method on two-dimensional passive 
advection test problems. 

In discussing unstructured grids, it is necessary to define the nomencl atures. To reitérate, the basic 
control volume element in our structured-unstructured computational domain is a truncated triangular 
prism. Each prism is bounded by five faces. For advection across each face, it is convenient to define a 
local coordínate system with its origin located at the center of the face. Each face separates the left- 
hand-side (LHS) from the right-hand-side (RHS) such that the flow from the LHS cell to the RHS cell is 
considered positive. For simplicity, the advected variable, hereafter denoted as \|/, is placed at the cell 
centroid, while the velocity vector is defined on the cell face at the origin of the local coordínate system. 
Figure 7 shows the basic arrangement of the variables on a two-dimensional grid. 



Figure 7. Arrangement of the variables on a 2-dimensional unstructured-grid. The generic 
advected variable, \|f, is placed at cell centroids. The subscripts L and R designates 
their left-hand-side and right-hand-side cell placements. The velocity is decomposed 
into normal and parallel components (not to scale). The flux F is co-located with V n 
at the face centers, and points only in the normal direction. 
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In its explicit form MPDATA adapts naturally to the above construct. As posed by Smolarkiewicz, the 
algorithm can be generalized to the following steps: 

1) At each cell face the low-order flux is found in conservative form using the standard first- 
order-accurate “upwind” scheme; 

2) The advected variable is integrated using the low-order flux; 

3) At each cell face, the low-order scheme is expanded in a Taylor series and the truncation 
error in the flux is explicitly identified; 

4) The error term is cast in the form of error velocity, V e ; 

5) The correction velocity, V c (= -V e ), is optionally limited to preserve monotonicity of the 
advected variable (Smolarkiewicz and Gradowski, 1990); 

6) Replacing V with V c , steps (1) through (5) are repeated a chosen number of times (=N C ) to 
achieve greater accuracy. 

Although Smolarkiewicz derived V c for a rectangular grid, it can be generalized to a grid with arbitrary 
control volume shape as long as the bounding faces are fíat. Consider the generic advection equation, 

^+V-(w)=0. (10) 

ot 

To compute the change in \|/ from time t = to to to+At, it is necessary to intégrate the flux V F / V / through 
'each face j during the period At: 

J 't= In+Aí — 

i ¡ dt a» 

Here a y = a ] é n is the area vector of face j, where é n denotes the unit vector normal to the face and 
pointing from left to right. For this integral, 2nd-order accuracy in space is achieved automatically by 
placing F = V y at the center of each face (in practice, at the faces is obtained by interpolation). 

Similarly, to ensure 2nd-order accuracy in time, ^V. should be evaluated at t = to+At/2. Assuming that 
a leapfrog algorithm is used (i.e., V,. is defined at t = to+At/2), we need only expand 'F. in a Taylor’s 
series as 

(^F At 

' ¥ i=' i 'j + ^r^- + °W’ ( 12 ) 

ot 2 

where the superscript 0 denotes an evaluation at t = to- Substituting (3.10) for V P / in (3.9) and 
performing the time integral, Equation 3.11 becomes 

<i3) 


Now substituting (3.10) for —- in (3.13) we have 

dt 
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(14) 


A'P; = 1%° - fv, • Vr, + (v ■ v, ]y }■ V, 

If we further let v n = V. • é n denote the component of the velocity normal to the face (and thus 
V. = v n e n ), then the first order upwind flux term is given by 

A^T*" =4v. +|v,||«»K +(v. -||v,|í> 8 }'a,A/ (15) 

Now letting S L denote the vector pointing from the cell centroid on the left to the face, and S R the 
vector pointing from the face to the cell centroid on the right, we can write 

v F i = 4 , ; -^-V^., (16a) 


and 


Similarly, we will let 


and 


+ SfW l 

W L = 'V r 'V L , 

9v, = 'v,-'r, 


(16b) 

(17a) 

(17b) 


Now we can rewrite (13) as: 


AT-- (18) 

The correction term is the difference between (14) and (18). After some algebraic manipulation, this 
correction term can be written as 


A'F, = 





rJ 


~ V n 


- V*F, /- - \ 
v .-=i-+(v-v.) 

1 VJ/ ' 1 > 


A t_ 
~2 


VajAl, 


(19) 


_ w + u/ 

where 'F = —-- L . 

2 

The correction velocity is now: 


( 'v aZ 1l 

4',+>P l 


-v. 


- VV, /- - X 


A t 

T 


v, 


v, 


( 20 ) 




(a) Computational mesh used for the 
rotating cone advection test. 


(b) Initial Conditions. 


Figure 8. Rotating cone advection test. 


In actual implementation, the correction term becomes 


A^ c = 


vñ+'F.y 


- V, 


V,. •V'F,. - - 

7 +V-V, 


A t 
~2 


^L,R a j^ > 


( 21 ) 


where V F is the valué of the advected quantity following the first-order upwind advection step. 

We limit the number of correction steps to N c =l, since additional correction steps are not cost effective. 
The effect of additional correction steps is to bring the solution closer to 2nd-order accuracy, which is 
very nearly achieved with just one correction step; the additional accuracy attainable is limited. 

To compare with known results (Smolarkiewicz, 1984, Smolarkiewicz and Clark, 1986, and 
Smolarkiewicz and Grabowski, 1990), we performed the rotating cone test on triangular grids. Figure 8 
shows the mesh and the initial placement of the passive scalar. Figure 9 shows the Solutions after six 
revolutions of the cone. Figure 9 shows the results of the Ist-order "upwind" scheme (N c =0). As 
shown, the original cone shape is completely smeared out. In Figure 9, we show the case with one 
correction step. This case has a background valué of zero, which when coupled with the positive 
definiteness aspect of the algorithm, generates extra numerical diffusion which automatically ensures 
monotonicity (Smolarkiewicz and Grabowski, 1990). Nevertheless, the cone shape has been 
substantially preserved. In general, all of our results are consistent with previous tests (Smolarkiewicz, 
1984, Smolarkiewicz and Clark, 1986, and Smolarkiewicz and Grabowski, 1990). 

The claimed near second-order-accuracy of our scheme is demonstrated by testing how the accumulated 
error in the computation scales with grid size. Three measurements of error are used for the rotating 
cone test: the root-mean-square error over the whole mesh, the root-mean-square error over an area 
extended by the original base of the cone, and deviation in the height of the cone. The mesh sizes range 








Figure 9. Rotating cone advection test after five rotations. 


from one and a half to one half of those shown in Figure 8, i.e., a factor of three variation. The results 
are shown in Figure 10. We note that except for the left-most data point, all three measurements of error 
scale approximately as A" (the dotted lines indícate the slopes of both A and A scaling for comparison). 
The left-most data points were generated from runs with the largest grid size. By observation, we find 
that the excess diffusion in these cases had reached the boundary, thus contaminating the runs results. 

3.2.2 Implicit Acoustic Filtering and Diffusion. 


One of the major goals of the OMEGA model development is to achieve a simulation capability that 
transcends scales through the use of variable resolution. The accurate portrayal of aerosols in the 
planetary boundary layer will require fine resolution in the structured vertical direction and selective 
locations of fine horizontal resolution. Accurate, stable numerical Solutions require that the time step of 
integration be determined by the finest grid spacing. This problem is magnified when using the fully- 
compressible set of nonhydrostatic equations, since fast-moving sound waves are inherent in the 
solution. Indeed, in the structured direction, the CFL time-step can become too restrictive. To 
overeóme this restriction, the growth and propagation of acoustic waves are controlled through the 
application of vertical and horizontal filters that act on the vertical component of the momentum. 

The momentum equation is integrated using a semi-implicit, fractional time step scheme. After explicit 
calculation of the advection and source terms, an implicit calculation is performed in which a correction 
term is added to the vertical component of momentum. This term, which has the effect of slowing the 
growth and speed of sound waves, is described below. The next pardal step consists of the application 
of a horizontal filter, which is described in the following section. Like the implicit correction, this term 
is based on and added to the vertical momentum. The implicit calculation of the eddy-diffusion term 
completes the momentum advance to time step m+1. 

The implicit calculation solves the following equation for w": 


w" = w" +C 5 V S 


di 


w n - w n 


dt 2 , 


( 22 ) 
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where w n is the vertical component of momentum at the beginning of the nth timestep and w" is the 
valué after the addition of the (explicitly calculated) terms representing advection, the Coriolis forcé, the 
pressure gradient contribution, and the buoyancy effect. Here, v s is the sound speed, and c s is an 
empirically chosen coefficient. A stability analysis of the above equation shows that c s must be greater 
than 1.0. The default valué currently used in OMEGA is 2.5. 

Eddy-diffusion is also treated implicitly in the radial direction. The ffactional step approach is necessary 
because the acoustic filter described above effectively increases the inertia at short wavelengths, i.e., 
Affective ~ p(l+k 2 ) where k is the wave number in the direction of propagation. This can reduce the action 
of the eddy-diffusion term unless the operations are split into sepárate fractional steps. 

Again, this implementation of the dispersión operator removes the severe time step restriction that can 
result ffom diffusion flux along the fine scale radial direction. A sepárate explicit treatment is used for 
momentum diffusion in the horizontal direction. 

3.2.3 Horizontal Acoustinc Filtering. 

The horizontal filter is based on the horizontal components of the laplacian of the vertical momentum. 
The coefficient of the horizontal filter is based on vertical height; no horizontal filtering is done at either 
the top boundary or the surface. The máximum valué of the coefficient is self-adjusting, initially having 
a negligible valué and increasing by small increments based on the growth of the laplacian of w. The 
Laplacian of the vertical component of the momentum is calculated using Green’s formula: 



Figure 10. Error scaling with grid resolution. Three measurements of deviation from 
the original cone are used: root-mean-square error over the whole mesh, 
root-mean-square error over an area extended by the original cone, and 
deviation in the height of the cone. 
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where A is the area enclosed by a curve. The discretized form becomes 


V 2 w 


4i 


dw 


A^á h, 


dL 


( 24 ) 


where A is now the area of the triangle, n, is the outward normal vector to the i‘ h edge, and dL ¡ is the 
length of the i ,h edge. The normal derivative is defíned as: 


dw = W,.-W c (2 5) 

ai,. d i 

where w¡ is the valué of w at the centroid of the triangle that shares the i' h edge of the target triangle, w c 
is the valué of w at the centroid of the target triangle, and d¡ is the distance between the two centroids. 

This approximation is first-order accurate only if the line segment connecting the centroids is 
perpendicular to the shared edge; i.e., this line segment coincides with the outward normal vector. To 
improve the accuracy of the approximation, we define a “ghost” triangle whose edges are the 
perpendicular bisectors of the line segments connecting the centroid of the target triangle with the 
centroids of its neighbors. The area and edge lengths of this ghost triangle are used in the formula for 
V 2 w. Figure 11 shows the ghost triangle constructed around the centroid of the inner cell. 


3.2.4 Horizontal Diffusion. 

The horizontal diffusion term in the momentum equation can be expanded as: 

V-(k¿Vm)=VK¿ ■Wü + K d 'V 2 ü (26) 

where K d is the horizontal diffusion coefficient, given by: 

K¿ = 0.36DA 2 . (27) 

Here D is the deformation rate, and A is the horizontal grid spacing. Assuming local deformation 
gradients are small, this can be written as: 

V-(K d V«)= K¿V 2 w . (28) 

We define orthonormal vectors 3c and y lying in a plañe tangent to the surface of the earth at the 
centroid of the triangle in which the deformation rate is sought (Figure 12). Projecting the velocity u 
onto this local coordínate system gives u and v, the horizontal components of velocity. By assuming 
the velocity to be locally linear in x and y , it is possible to express u and v as a u x + b„y + c u and 
a v x+b v y +c v , respectively. For the horizontal deformation rate D we will use the approximation 

Using the valúes of u and v at the vértices of the triangle, it is then possible to solve for b„ and a v , 
which gives the local deformation rate as D = b u +a v . In this way, actual horizontal deformation is 
obtained, as opposed to deformation in the plañe of the triangle. Since the horizontal grid spacing is 
squared in the expression for the deformation coefficient, we use the area of the triangle 




Figure 11. The ghost triangle used in the calculation of the Laplacian. 

D = (29) 

ay dx 

The Laplacian of each component of the velocity is calculated using the discretization of Green’s 
formula based on the “ghost triangles” described in the previous section. 

3.2.5 Model Boundary Conditions. 

3.2.5.1 Surface Boundary Conditions . The surface bnoundary conditions applied in the OMEGA 
* model are based on a diagnostic layer treatment, which is Dirichlet in its implementation but has a 

Neumann flavor, since new surface valúes are calculated based on flux estimations. Surface boundary 
conditions are discussed more fully in the section on planetary boundary layer physics. 

3.2.5.2 Lateral Boundary Conditions. In cloud scale models, in which the time scale is small enough 
to assume that the environment outside the model domain is unchanged during the simulation, inflow 
boundaries are usually prescribed from environmental conditions. In mesoscale and larger models this 
assumption can no longer be made. The environment outside the computational domain has to be 



Figure 12. A local coordínate system defines the horizontal plañe 
passing through the centroid of each triangle. 
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derived from larger scale forecast tools such as the LFM (Limited Fine Mesh model), NGM (Nested 
Grid Model), ETA model, or MRF (Médium Range Forecast model) run at the National Centers for 
Environmental Prediction (NCEP). Boundary conditions at user-specified time intervals are calculated 
based on one of these forecasts, and linear interpolation is used to determine boundary valúes at 
intermedíate times. 

The OMOEGA model utilizes “ghost cells” at all boundaries. This single layer of cells holds predicted 
valúes of all prognostic variables obtained from the chosen forecast data. Boundary faces are treated in 
the same way as interior faces throughout the calculation, with one notable exception described below. 
At the end of each time step, a “nudging” scheme is applied to the new valúes of all variables in cells 
immediately inside the boundary to the forecast valúes in the ghost cells across the boundary. 

It is worth noting that linear interpolation of boundary conditions between two forecast times can fail to 
accurately reflect the diumal cycle of heating and cooling that affects temperatures and winds near the 
surface. This is particularly true if data is available only at twelve-hour intervals. If, for example, data 
is available for a certain lócale at 7am and 7pm, the mínimum and máximum temperatures obtained in a 
twenty-four hour period will not appear in the boundary data. This will lead to false temperature and 
pressure gradients across lateral boundary edges near the surface, which we ignore by assuming zero 
pressure gradients across lateral boundaries for the purpose of calculating momentum. Thus we avoid 
the inducement of sea-breeze-like or off-shore-breeze-like winds along lateral boundaries. This is not an 
ideal solution; nudging of the prognostic variables along lateral boundaries still corrupts temperatures 
and pressures during certain portions of the diumal cycle, however, the error is not cumulative, ñor does 
the affect extend far into the domain. 

3.2.S.3 Toy Boundary Conditions . At the top of the computational domain, initial valúes of density 
and potential temperature are calculated assuming a hydrostatic balance, and these valúes remain fixed 
throughout the simulation. A homogeneous Neumann boundary condition is applied to the horizontal 
components of velocity, while the vertical velocity is assumed to be zero. Thus the upper boundary is a 
rigid, free-slip surface. To minimize reflection off of this boundary, a diffusive “sponge” is applied to 
the vertical component of the velocity. The coefficient of the sponge is exponential, and approaches 
zero rapidly below the top few layers of the grid. 

3.2.6 Time Splitting. 

The OMEGA user interface and grid generator allow the user to specify selected regions in which the 
resolution can be as fine as one kilometer or less. The cell size in the outer portion of the domain may 
be more than two orders of magnitude greater than this. For this reason, it is highly desirable to utilize a 
time-splitting scheme in which the hydrodynamic solver is sub-cycled over small cells, thus preventing 
the time step over the entire domain from being limited to that required by the smallest cell. 

In OMEGA, we use the term “inner mask” to refer to a collection of cells over which the solution will 
be sub-cycled. This terminology distinguishes such a collection of cells from a “subdomain,” which is 
user-specified as part of the grid definition process. An inner mask is selected by the model, which 
calculates the Courant-limited time step for each cell based on sound speed, and uses this as the criterion 
for mask definition. Cells in an inner mask need not be contiguous. We impose only one restriction: 
any cell having two neighbors in an inner mask is added to the inner mask. For stability, no more than 
two time steps are taken in a cell before the neighboring cell is updated; however, masks can be nested, 
so that the time-step in the outermost portion of the domain can be many times that of the smallest cells. 

The difficulty in implementing a time-splitting scheme lies in the treatment of the mask boundaries. In 
the OMEGA scheme, a face shared by cells in different masks belongs to the innermost mask. Thus 
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advection across this face takes place at the smaller timestep associated with the inner mask. Advected 
quantities are updated in the cell belonging to the inner mask, while the change in the quantity on the 
cell in the outer mask is stored. Advected quantities in the outer cell are updated when advection across 
faces belonging to the outer mask occurs. In this way, the accuracy of the advection scheme is formally 
preserved for both cells, each at their respective time step. This process is aided by the application, at 
each sub-step, of the horizontal filter to a larger mask consisting of the inner mask and its bounding 
cells. This filter, which dissipates horizontally traveling acoustic waves, helps to maintain a smooth 
solution across mask boundaries. In numerous test cases involving as many as five nested masks, this 
time-splitting scheme has proven to be stable and robust. 

3.2.7 Model Validation. 

In the course of developing the OMEGA model, we have assembled a set of test cases that are used to 
validate various aspects of the model. We continuously rerun these cases each time a significant change 
is made. In addition to a number of analytic cases which test specific aspects of the solution, we rely on 
a collection of real cases for which we have done extensive analysis of the available meteorological and 
diffusion data. 

The simplest analytic test cases consist of puré advection of a passive scalar in either a simple rotational 
flow field or a more complex flow field consisting of counter-rotating vórtices. Figure 13 shows the 
result of a passive scalar in a steady State flow field composed of six counter rotating vórtices. A second 
order accurate solution of the OMEGA model solution after 1000 time steps (right panel) agrees well 
with the corresponding analytical solution(center). 

The time splitting technique used in the OMEGA model is tested with a rotating cone test case discussed 
earlier in this section. The domain of integration contains three nested "masks" or (Figure 14) over 
which the hydrodynamics are subcycled to permit a larger timestep on the outer portion of the domain 
and smaller time steps in the inner masks. The cone initially located over the Southwest región of the 
domain (Figure 15) rotated clockwise during the model integration and passed over the high resolution 



Figure 13. Advection test consisting of a passive scalar in a steady-state flow 
field composed of six counter-rotating vórtices (left). The analytic 
solution at a given scaled time (from Staniforth) is shown (center); 
along with the second order solution at the same scaled time (1000 
timesteps) for the OMEGA advection solver (right). 













Figure 14. Coarse and two finer subdomains for testing time splitting technique. 


región and moved over the southeast región (Figure 15). A comparison of these two figures shows that 
the cone has not been distorted by the use of time splitting technique. A detail discussion of these runs 
and validation runs against meteorological and diffusion data sets are given in Volume 2 of this final 
report. 



(a) Initial position of the cone. (b) Final position of the cone. 

Figure 15. Initial and final positions of the cone. 
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Section 4 


Model Physics 


This section discusses the five most important processes to atmospheric circulation: turbulence, air 
surface interactions, cloud microphysics, convective parameterization, and radiation transport. 

Turbulence affects, among other things, the rate at which surface moisture and energy enters the 
atmosphere. Cloud microphysics and convective parameterization are important in the vertical 
distribution of this moisture and energy and in the formation of precipitation. And finally, radiation 
transport is important in determining the solar heating of the Earth’s surface and the atmosphere. 

4.1 Turbulence and the Planetary Boundary Layer. 

Much has been written in recent years conceming modeling and parameterizing atmospheric turbulence. 
However, little has been resolved in terms of the development of improved parameterizations or the 
direction and future progress of turbulence parameterization schemes as they relate to numerical weather 
prediction. The parameterization of turbulence in OMEGA (i.e., the forcing terms, F, in Equations 3,4, 
and 5) are divided into two parts: horizontal and vertical. The horizontal diffusion is a function of the 
deformation of momentum (Section 3.2.2), while the vertical diffusion is implemented using a multi- 
level planetary boundary layer (PBL) model. 

The effects of the PBL can be incorporated into a mesoscale model in two ways. One way is to 
parameterize the entire PBL as one layer. This involves identifying and relating unresolved processes in 
the PBL with resolvable ones. The complexity of this single-layer PBL parameterization lies in the 
variety and interdependence of atmospheric processes acting on different scales. The second approach is 
to inelude several computational levels in the PBL in order to resolve the boundary layer structure 
effectively and explicitly. Such multi-level PBL formulations require turbulent fluxes of momentum, 
heat, and moisture at several levels within the PBL. Thus they require some type of closure scheme to 
relate turbulent fluxes to mean quantities. The atmospheric PBL in OMEGA is treated separately as the 
viscous sublayer, the surface layer, and the transition layer using a turbulent kinetic energy (TKE) 
closure scheme. 

4.1.1 Viscous Sublayer. 

The viscous sublayer is defined as the level near the ground (z < z 0 ) where the transfer of the dependent 
variables by molecular motion becomes important. Temperature and specific humidity in this layer are 
related to their ground surface valúes (Deardorff, 1974): 

e z0 =e c + 0.0962 (0. /k){u. z 0 /v) 045 , (30) 

and q z0 = q G +0.0962 {q./k){u . z 0 /v) 045 . (31) 

By analogy for other gaseous and atmospheric materials (%) 

Z* =Xc +0.0962Ctr. /*)(«. z 0 /v) 0 45 • (32) 

In these expressions, v is the kinematic viscosity of air, k , von Karman's constant, 14 , the friction 
velocity, z 0 , the surface roughness length, 6*, the subgrid scale temperature flux, and q t , the subgrid 
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scale moisture flux. In this layer, it is assumed that u = v = w = 0, and the variation of pressure is 
ignored. 

4.1.2 Surface Layer. 

There is no precise definition of this layer, (which may contain more than one model layer) but it is that 
part of the planetary boundary layer immediately above the viscous sublayer where small-scale 
turbulence dominates the transfer of energy, moisture, and momentum, and the vertical variation of the 
vertical fluxes is roughly less than 10%. For this reason it is also referred to as the constant flux layer. 
The depth of this layer usually varíes from about lOm to lOOm. It has been of considerable interest to 
micrometeorologists to find a suitable theoretical or semi-empirical framework for a quantitative 
description of the mean and turbulence structure of this layer. One of the most common formulations 
used in mesoscale models was reported by Beljaars and Holtslag (1991) and in which 


u.=kVj [ln(z ¡/ z 0 )-Y M fo /L) + Ym (z 0 1 L )\ - ( 33 ) 

6. = k (6 X -6 s )¡ [ln(z,/z 0 )— Yh fc l L ) + Wh (z 0 ' L )]> ( 34 ) 

q* =k{q x -^)/[ln(z 1 /z 0 )-V' w (z,/¿)+^//(z 0 /L )]’ ( 35 ) 

and X*=k{x x -Xs)l\^(zxlza)-YH{z\l L )+Y H {z Q IL)], (36) 


where the symbols V,, 6 X , q x , and X\ are wind speed, potential temperature, specific humidity, and 
other gaseous concentrations at the first model level in the atmosphere, 0 S , q s , and x s are the 
corresponding surface valúes, u, is the friction velocity, 9, is the temperature scale, q. is the humidity 
scale, x * is the other gaseous scale, k is the von Karman constant (k = 0.4), is the height of the first 
model level, and L is the Monin Obukhov length, Ym and Yh are dimensionless stability functions of 
height divided by the Monin Obukhov length. 


These integral functions for unstable surface layer are given by 

Y m = 21n[(l + x)/2] + ln[(l + jc 2 )/ 2 J — 2 tan -1 (jc) -H / 2 , 
Y h =21n[(l + x 2 )/2], 

where 


jc = (l -16 z!L) v \ 

In stable surface layer the following relationships are used 


Ym =a- + b\ 
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where a=l, b=0667, c=5, d=0.35, and z/L is a non-dimensional stability parameter. The Monin- 
Obukhov stability length, L, is defined as 

0i¿ 

kg0.’ 

where 6 is the average air temperature near the surface, and g is the gravity. 


( 42 ) 


4.1.3 Transition Layer. 

The turbulent región extending above the surface layer up to the top of the planetary boundary layer (i.e. 
90% of its depth) is called the transition layer. Its depth is variable since it depends on the strength of 
turbulence generated at the ground surface. The PBL height often can be considered to be the depth of 
turbulent activity. Therefore, defining the top of the PBL is not always simple. For example, on a sunny 
day with light winds thermal convection may push this layer up to l-2km above the surface. At night 
with weak winds and little cloud cover it may contract to less than lOOm. One of the most important 
characteristics of this layer is that the effects of roughness in generating forced convection decrease more 
rapidly with height than effects of heating on thermal convection. Therefore, this layer tends to be well 
mixed and is generally characterized by free convection involving rather large eddy sizes associated with 
thermals and heat plumes. Turbulence associated with these large eddies can be parameterized in terms 
of the mean flow. This leads to the closure problem that is addressed in this section. In OMEGA a 
turbulent kinetic energy closure scheme is used. 

4.1.3.1 Turbulent Kinetic Enersv (TKE) Closure. The turbulent kinetic energy (TKE) closure scheme 
is an improvement on the relatively simple first-order closure scheme. TKE closure accounts for more 
of the physics of the atmosphere in the formulation of the eddy diffusivity, K. TKE closure begins with 
the same fundamental assumptions of first-order closure, but in the TKE method the K coefficient is 
determined from more physically realistic relationships given by prognostic equation for turbulent 
kinetic energy. This type of closure is more economical than higher-order closure schemes. For 
example, second order closure requires the integration of at least 15 partial differential equations when 
moisture is included. Second-order closure also requires initial and boundary conditions for terms such 
as the Reynolds Stress that are not operationally available; this, coupled with the computational expense 
renders second order closure inappropriate for an operational model such as OMEGA. 


In contrast with second order closure, TKE closure requires only one additional equation, yet it provides 
a more physically realistic solution to the closure problem than does first-order closure. The vertical 
turbulent exchange is parameterized using the “2.5 level” closure technique developed by Mellor and 
Yamada (1974). The level 2.5 model has been selected from the hierarchy of available TKE closure 
schemes because it is a good compromise between accuracy and efficiency, and because it has been 
applied to many atmospheric modeling studies. The level 2.5 turbulence closure is based on the 
following prognostic equation for the turbulent kinetic energy e: 


¿(pe) 

dt 


= -V.(pev) + V r (K e V r pe) 


+ 


-u w 


d u 
dz 


■dv 
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vV—+ 7w’0’ 

dz e 


-e 


(43) 


The first two terms on the right hand side represent advection and turbulent transpon of TKE (V r 
represents the radial gradient), the third term in the bracket represents both shear and buoyancy 
production of TKE, and the fourth term dissipation of turbulent energy. In this closure scheme, several 
terms on the right hand side must be parameterized. Following Deardorff (1980), the eddy fluxes are 
parameterized as 
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(44) 


W ’0’ = -K h ^p-, 

h dz 

(45) 

and 


(46) 


The parameterization of the final term, energy dissipation e , is achieved diagnostically using 
Kolmogorov’s hypothesis as 

e = c £ e m /¿, (47) 

where c £ is a constant assumed to be 0.17. 


In this closure scheme, the prognostic equations of all other second-order moments are reduced to 
algebraic equations. It therefore allows the expression of eddy diffusivities after considerable algebraic 
reduction in the following simple form: 



<s 

N 

II 

(48) 


K h =t(2e) m S„ , 

(49) 

and 

K,=e(2e)' n S„ 

(50) 

where non-dimensional eddy diffusivities S m and S h are obtained from 



^ m S fj (l 3 A 2 B 2 G h 12 A l A 2 G H ) — A 2 

(51) 



(52) 

(53) 


= _^8_M v 

2e 6 dz 


( 54 ) 


In the above equations, empirical constants are assigned valúes after Mellor and Yamada (1982), 
(A 1 ,A 2 ,B 1 ,B 2 ,C l ,S e ) = (0.92,0.74,16.6,10.1,0.08,0.20). It is also assumed that G H < 0.033 and 
Ga, < 0.825 - 25 G h to reduce the singularities of the level 2.5 closure. 


The turbulent length scale l is still unknown and needs to be parameterized in terms of the known 
variables. After Mellor and Yamada (1982) 
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( 55 ) 


t _ k{z + z 0 ) 
l+k{z + z 0 )/e„ 

where £ n the master (or primary) length scale defined as 

H 

jjezdz 



(56) 


with a, a constant coefficient of the order of 0.1. In order to smoothly change the mixing length from 
the neutral to the slightly stable condition, the stable mixing length is modified with an upper limit as 
proposed by Andre et al. (1978) such that 


l < 0.75 


le 


8 

Ld dz 


1/2 


(57) 


4.2 Air Surface Interations. 

OMEGA requires as boundary conditions the temperature and humidity at the roughness level. These 
must be calculated prognostically taking the interaction of the atmosphere and the land surface into 
account. Numerical studies have demonstrated that contrasts in land surface characteristics (e.g., soil 
texture and moisture, vegetation type) can generate local circulations as strong as sea breezes (McCorcle, 
1988; Chang and Wetzel, 1990). In consequence, landscape variations strongly influence atmospheric 
dispersión pattems. 

The land surface module in OMEGA is based on the scheme proposed recently by Noilhan and Plantón 
(1989). Following their approach, only one energy balance is considered for a whole soil-vegetation 
system. A single surface skin temperature T s , and surface air humidity q s are computed for both the 
canopy and the ground. The calculation of two distinct temperatures for the soil surface and the 
vegetation would significantly increase the complexity of the scheme since it requires an evaluation of 
radiative and turbulent transfers between the soil, the plant, and the air inside the canopy. 

The implemented scheme allows for two layers in soil and a single vegetation canopy. It computes the 
evolution of four variables: the surface temperature of the soil-vegetation médium T s , the soil water 
contents rj l and r¡ 2 in the upper and lower soil layers, and the amount of liquid water retained on the 
foliage W d . 

The parameters used to described the main physical processes in a soil-vegetation system can be divided 
into two classes: primary parameters to be specified by spatial distribution, and secondary parameters 
whose valúes can be associated with the valúes of the primary parameters. The land-use in the modeling 
domain is classified into several categories characterized by two numerical indices which are the primary 
parameters: type of soil texture and type of land. The secondaiy parameters describing a variety of 
physical characteristics of soil are listed in Table 1. 
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4.2.1 Surface Temperatura. 

The most indicative variable of surface variation is ground temperature T s . OMEGA calculates ground 
temperature using the force-restore rate equation of Deardorff (1978) 

~rr 

-£ = C t G 0+ Cí(T Js -T,), (58) 

where Q is angular velocity of earth (fi = 7.27229xl(T 5 s~' ), and G 0 is the heat storage rate in the soil- 
vegetation médium, which is equal to the sum of all the atmospheric fluxes at the surface 

«o = K + »«*, -(H 0 +L t E 0 ), (59) 

where R 0 is the net radiation at the surface (see next section for detail), H anthrop is the anthropogenic 
heat flux, H 0 and L C E 0 are the sensible and latent heat fluxes from the atmosphere. Surface turbulent 
fluxes are parameterized using Monin-Obukhov surface-similarity theory. The sensible heat flux to the 
atmosphere, H 0 , the latent heat flux from evaporation or dew formation, L C E 0 , are calculated using the 

turbulent subgrid scale friction velocity, «*, temperature scale, 0»,and humidity scale, q t , 

H 0 = pC p m, 8, (60) 

and L C E 0 = pL c u.q* . (61) 

The forcing terms in the force-restore equation may be completed by specification of the anthropogenic 
heat flux H an!hrop if urban or industrial areas are considered. 


The first term on the right-hand side of Equation (58) represents the diumal forcing of T s by the heat 
storage rate in the soil-vegetation médium, G 0 , and the second one tends to restore T s to the mean soil 
temperature T ds . 


The thermal capacity C T combines the capacities of both vegetation C , and soil C g weighted by 
vegetation cover (veg) and bare soil cover (1 -veg) (Noilhan and Plantón, 1989) 

1 _ veg | 1 -veg 

C T c c 

T ^veg '■'g 


(62) 


C veg has a valué (lO 3 AT l m 2 J 1 ) related to the weak heat storage and transfer ability of vegetation. The 
thermal capacity of soil 



(63) 


is given by its thermal inertia 


h=4Ps c s^s =—7 

V T i 


(64) 
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where p s is the density of soil, c s is the specific heat of soil, and X s is coefficient of thermal soil 
conductivity. T ds is the approximately constant mean deep soil temperature at the diumal temperature 
cycle penetrative soil depth 

d - = í¥^ (< 

\PsCs 

where t, is the period of 1 day. 

The soil properties in Equation (63) depend upon both the soil texture and the soil moisture. This 
dependence is deduced from the analytical expression given by (Noilhan and Plantón, 1989) 




b¡ 21ogl0 


where C , ij s , and b are estimated for each soil texture from table (4.1). 


4.2.2 Surface Fluxes of Radiation. 

At the surface, the net shortwave solar radiation R s (radiative contribution) is obtained by using two 
empirical functions. The first is a Rayleigh transmission function that ineludes the effeets of the 
molecular scattering and absorption by permanent gases such as oxygen, ozone, and carbón dioxide 
(Atwater and Brown, 1974) such that 


G = 0.485+ 0.515 1.041-0.16 


0.000949/7 + 0.051 
cosZ 


where p is pressure in mb. A second function accounts for the absorptivity of water vapor (McDonald, 
1960) 


r(z) 

= 0.077 , 

_cosZ_ 


where r is the optical path length of water vapor above the layer z given by 

top 

r(z)= ¡pqdz . (< 

The net shortwave solar radiation at the surface is 

R s ={ ~ X^I cosz + «2 - «wXl - ) eos Z>-a 2 /a l 
[ 0 cos Z<-c 2 /a, 

where A g is the albedo, the empirical coefficients are a¡ = 990 Wm~ 2 and a 2 = -30 Wm ~ 2 , and Z is the 
zenith angle 

cosZ = cos^ cos S cos H + sin y/ sin¿> . (' 

íí^is the latitude, 8 the solar declination, which is a function of Julián day, and H the solar hour angle 
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( 72 ) 


u _ n{t-\2) 

12 

where t is local time in hours. 

In the above equation, the cloud effects on the incoming radiation is also parameterized in terms of the 
effective cloud cover N E following Kasten and Czeplak (1980) 

N e ={N t +N l )/2, (73) 

where N T is the total cloud cover, and N L is the cover due to low clouds. The remaining empirical 
coefficients are ¿>, =0.75 and b 2 - 3.4. 

The net longwave radiation has two components: (1) incoming radiation from the atmosphere 
(atmospheric radiation), and (2) outgoing radiation from the surface (terrestrial radiation) 

R^R LÍ -R Lt . (74) 

The incoming longwave radiation is expressed by 

R Ll =c ¡ T‘+c 2 N £ , (75) 

where c¡ = 5.3 brt0~ 13 Wm~ 2 K ~*, c 2 = 60 Wm~ 2 , and the reference air temperature T a at z=50 m is 
recommended by Van Ulden and Holtslag (1985). The outgoing longwave radiation from the surface 
arises from the Stefan-Boltzmann law 

R l7 =£ s cT s \ (76) 

where £ g is the emissivity for the ground surface specified in terms of land-use categories, <7 is the 
Stefan-Boltzman’s constant (5.67xl0' 8 Wm' 2 K 4 ). 

Finally, the net radiation flux at ground level is defined as 

R 0 = R s +R l . (77) 

The surface energy budget, (Equation 59), is updated each timestep to provide an updated ground 
temperature. The ground temperature is then related to the potential temperature in the lowest cell by the 
following equation: 


dsfc Tg 


P 

. P ) 


where p re /is the reference pressure, 1000 mb. 


(78) 


4.2.3 Soil Hydrology. 

The soil wetness influences the energy balance of the surface by changing the surface albedo by affecting 
the soil heat flux through the thermal diffusivity, and by determining the partitioning of the heat fluxes 
into latent and sensible heat flux. The 2-layer soil hydrology model developed by Mahrt and Pan (1984) 
is implemented to represent the soil moisture. 
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The transport of water in the soil can be described by the prognostic equation 

dr¡ _ 1 dW s 
dt p w dz 


■ + A(z), 


( 79 ) 


where tj is the volumetric soil moisture content, p w is the density of liquid water, W s is the moisture 
flux within the soil. The term A(z) describes extraction of the soil water by roots due to the 
transpiration E tr from vegetation canopy. 


The soil moisture flux is expressed as 


or in altérnate form as 


= Kn 


¿QF + z) 

dz 


W s = D n p w ~ + K J¡Pw , 


( 80 ) 


(81) 


where K n is the hydraulic conductivity, 'F is the moisture potential, which represents the potential 
energy needed to extract water against capillary and adhesive forces in the soil, and the diffusivity, D n , is 
defined as 

d¥ 


D ’- K *>n 


(82) 


The parameters K T¡ ,D TJ , and *F are related to tj through a set of simple relationships presented by Clapp 
and Homberg (1978) and applied in McCumber and Pielke (1981) 


= K m 
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(83) 


(84) 


(85) 


The exponent b, saturated hydraulic conductivity K m , saturated moisture potential 'Fj, and soil moisture 
at field capacity (soil porosity) r¡ s are specified as function of the USDA soil textural classes determined 
by their percentage of clay, silt and sand (Table 1). 
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Table 1: Soil parameters of the USDA (United States Department of Agriculture) textural classes. 



Type of soil 

Vs 

Vwiilt 


Vs 


b 

P¡ c i 

A) 


1 

Sand 

0.395 

0.068 

0.135 

-0.121 

1.76X10" 1 

4.05 

1.463x10 6 

0.26 

3.222 

2 

Loamy sand 

0.410 

0.075 

0.150 

-0.090 

1.563X10" 1 

4.38 

1.404x10 6 

0.27 

3.057 

3 

Sandy loam 

0.435 

0.114 

0.195 

-0.218 

3.41xl0' 3 

4.90 

1.320x10 6 

0.28 

3.560 

4 

Silt loam 

0.485 

0.179 

0.255 

-0.786 

7.2x10‘ 6 

5.30 

1.271xl0 6 

0.29 

4.418 

5 

Loam 

0.451 

0.155 

0.240 

-0.478 

7.0xl0 ¿ 

5.39 

1.212xl0 6 

0.30 

4.111 

6 

Sandy clay loam 

0.420 

0.175 

0.255 

-0.299 

6.3x10' 6 

7.12 

1.175xl0 6 

0.28 

3.670 

7 

Silty clay loam 

0.477 

0.218 

0.322 

-0.356 

1.7xl0' 6 

7.75 

1.317xl0 6 

0.26 

3.593 

8 

Clay loam 

0.476 

0.250 

0.325 

-0.630 

2.5X10 6 

8.52 

1.225x10 6 

0.24 

3.995 

9 

Sandy clay 

0.426 

0.219 

0.310 

-0.153 

2.2x10"* 

10.4 

1.175xl0 6 

0.22 

3.058 

10 

Silty clay 

0.492 

0.283 

0.370 

-0.490 

1.0x10' 6 

10.4 

1.115xl0 6 

0.21 

3.729 

11 

Clay 

0.482 

0.286 

0.367 

-0.405 

1.3x10' 6 

11.4 

1.089x10 6 

0.20 

3.600 

12 

Silt 

0.485 

0.179 

0.535 

-0.786 

8.0x10' 6 

5.30 

0.836xl0 6 

0.14 

3.970 


Following Mahrt and Pan (1984) equation (4.50) is integrated over two slabs of the soil, from z=0 to 
z=-0.05m, and from z=-0.05m to z=-1.0m. This gives us two prognostic equations for the mean 
volumetric water content in the upper thin layer, 77 ,, and in the lower later, r¡ 2 


dt 


Zx =- 


D, 


dri 

dz 


P -E 




P VI 


p 

root Ir * 


( 86 ) 


*lH Z - 7 )=-(d ^ 
dt [Zl Zl) \ ' dz 


+ K ni -K, 2 -(l-f roo ,)E Ir . 


(87) 


The root extraction term A(z) was integrated under assumption that the upper soil layer contains 
fraction f roo , of root mass while the remaining ffaction of roots 1 - f wot is within the lower soil layer. 
Continuity at the soil surface requires that the flux W s of water in the soil at z = 0 equals of the sum of 
the flux P g of liquid water reaching the ground surface and the evaporation E g from the soil surface 


D ^ v - Eg 

k *~ p. 


( 88 ) 


The gradient of soil water content at the interfacial level, - z,, is evaluated as a finite difference 
assuming that the mean valué of 77 for each layer is represented for the midlevel of the layer 


f 

v 



)< 


= 2 D„ 


z 2 


(89) 


This term is taken equal (except sign) for both soil layers to maintain continuity. The diffusivity D n¡ 
and the hydraulic conductivity K n¡ are evaluated using the largest valué of 77 , and r ¡ 2 . The gradient 
term in Equation 89 is calculated between the surface and the midlevel of the upper layer using the 
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surface valué of soil water content Tj x and JJ 2 . It allows us to calcúlate diagnostically the near-surface 
soil water content r¡ x from Equation 86. 


Thermal properties of soil and its albedo vary depending on the amount of water present. The soil heat 
capacity, p s c s , is calculated as 

Ps c s={ l -n s )p i c ¡ +VP w c w , (9C 

where dry volumetric soil heat capacity, p i c ¡ , is tabulated for each soil class in Table 1, and p w c w are 

corresponding valúes for liquid water. The thermal conductivity of soil is given by McCumber and 
Pielke(1981) 


X s = max (419 ex p(-P / -2.7l),0.172), 


where P f = log(-'F/100). 


The valué of the soil albedo is determined taking into account the effect of soil wetness 

A g =A. n +A z , 


A = A,, (l - A), A = min {fi ,0.5), 


and the effect of solar zenith angle 


A z = 0.0l(exp(0.003286Z 1,5 )-1). 


Valúes of the basic soil albedo, Aq , are listed in Table 1, and the Halstedt parameter J3 g is defined by 

P = j0.5[l-cos(^7 0 / 77 /c )J T] 0 <Tj fc (95) 

g { 1 otherwise 

This expression takes into account that the surface evaporates at the potential rate when the soil moisture 
r¡ 0 exceeds the field capacity r¡ fc which corresponds to the hydraulic conductivity of 0.1 mm/day 

(Table 1). 


4.2.4 Water Intercepted by Canopy. 

The amount of liquid water, W d , intercepted by the foliage due to dew formation or rainfall is govemed 
by the prognostic equation (Deardorff, 1978) 

pm 

~r~ = vegP-E r + R r , (96) 

at 

where P is the precipitation rate at the top of the vegetation, E r is the direct evaporation from the 
fraction 8 of the foliage covered with a film of water when positive and the dew flux when negative. R r 



is the runoff of the interception reservoir which occurs when W d exceeds a máximum valué W dmax 
depending on the density of the canopy (Dickinson et al., 1986) 

W dtmi = KvegLAI, (97) 

where h w is the máximum storage of water, assumed to be equal 0.0002 m. The wet fraction of the 
canopy S is calculated following Deardorff (1978) 

forE r > 0 
[ 1 forE,< O' 

The flux of liquid water reaching the ground surface 

P g =(l -veg)P + R„ (99) 

ineludes the nonintercepted rainfall and the runoff from the canopy. 


4.2.5 Water Areas. 


The surface of water bodies such as lakes and oceans is represented in the OMEGA model rather simply. 
Because of the large heat capacity of water and strong surface mixing, the surface temperature of water 
bodies T sea is initially specified based climatological data and then assumed to be constant during the 
period of the mesoscale simulation. As in most mesoscale models, OMEGA assumes that the air just 
above the water surface is fully saturated by water vapor 

4o =4„(r-.) (100) 

The spatial variation of surface roughness length z Q is specified based on land-use categories over land 
surfaces, while it is determined using the Chamock’s relationship (Chamock, 1955) over water surfaces 


«o-«i 



( 101 ) 


Although the coefficient a x in the above equation has been found to be somewhat variable, a valué of 
0.018 is included in the model (Arya, 1988). Also, a mínimum valué of 0.0015cm is imposed on the z 0 
based on observations. 


4.3 Microphysics. 

Cloud microphysics is the general term used to represent the collection of processes that involve water 
and the different phases of water that occur in clouds. In nature, these processes inelude condensation / 
evaporation (vapor liquid), deposition/sublimation (vapor solid), freezing / melting (liquid <h> 
solid), interparticle collisions and collection, and sedimentation (fall out). Even though these processes 
appear conceptually straight forward, developing a complete mathematical formulation describing them 
is a very difficult task. Researchers subdivide cloud microphysics into several conceptual processes, and 
develop parameterizations for them guided by the actual physical processes they represent and with the 
relevant parameters derived from experimental and field observations. 

The OMEGA microphysics package is derived from that contained in the SAIC versión of the MASS 
model (originally designed by scientists at MESO, Inc.) and is a subset of the scheme described in Lin et 
al. (1983). These parameterizations fall under the category of bulk-water microphysics in which the 
production rates are functions of the total mass density of each water species. As in MASS, the water 
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Figure 16. Microphysical processes parameterized in OMEGA. 

species in OMEGA are divided into vapor, cloud droplets, ice crystals, rain and snow fields; however 
liquid and solid phases of water are not allowed to coexist. This implies that neither supercooled liquid 
ñor melting solids (such as wet snow) can exist in the model domain. Because of this simplification, the 
number of microphysical processes to be dealt with is much smaller than that in the cloud model TASS 
(Proctor 1987). The microphysical processes included in OMEGA (Figure 16) are: 

1) Condensation (or evaporatíon) of cloud water, 

2) Autoconversion of cloud water to rain, 

3) Condensation on (or evaporatíon of) raindrops, 

4) Accretion of cloud water by rain drops, 

5) Accretion of cloud water by rain drops, 

6) Fallout of rain, 

7) Generation of ice crystals from ice nuclei, 

8) Deposition (or sublimatíon) growth of ice crystals, 

9) Autoconversion of ice crystals to snow, 

10) Deposition on (or sublimatíon of) snow, 

11) Accretion of ice crystals by snow, 

12) Fallout of snow. 
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4.3.1 Computational Treatment. 

In OMEGA the time step used for integration depends on the local Courant number, which is a function 
of the velocity (advective or acoustic) in the cell and the dimensión of the cell. Henee the time step can 
range from sub-second to tens of seconds depending on the flow characteristics. The microphysical 
parameterizations, however, only need to be calculated using a time step of at most a few seconds. 

Henee the microphysical calculations are sub-eyeled in time if necessary. The microphysical time step is 
expected to be at roughly five to ten seconds. In the following discussion all valúes are assumed to be in 
SI units unless explicitly specified otherwise. In OMEGA the water substances are categorized into fíve 
species: vapor, cloud droplets, ice crystals, rain, and snow (hail and its associated physics could be 
included in the future). In the following discussion Q will represent the mixing ratio of a hydrometeor 
class. Suffixes V, C, I, R, and S will refer to water vapor, cloud droplets, ice crystals, rain and snow, 
respectively. 


The cloud droplets and ice crystals are assumed to be monodispersed and non-precipitating. Marshall 
and Palmer (1948) type relationships are used in describing rainfall droplets and snow crystals with their 
sizes. These are given as: 


and 


M¿*) = Wo* ex P ' 

V 

/ 

W 0 5 eX P ■ 

v 



( 102 ) 


where d R and ds are the particle diameters for rain and snow respectively. N 0 r and N 0 s are the y- 
intercepts of the distribution of rain and snow and henee represent the limiting number concentration as 
diameter tends to zero. A r and As represent the slope of the respective inverse-exponential distribution. 
The computation of the parameters and variables involved in the microphysical processes can be divided 
into several levels: 

• Parameters that are constants in space and time can be specified at the time of model 
compilation or at model run time. 

• Parameters that are functions of space but are constant in time can be precomputed as part of 
model initialization. 

• Parameters that are functions of space and time need to be evaluated in MICPHY (the 
subroutine in which the microphysical conversión rates are computed) for each time step. 
Some parameters depend only on a single physical variable which in tum is a function of 
space and time. For computational efficiency, look-up tables can be created in which these 
parameters are evaluated over a range of the dependent variable that will sufficiently bracket 
the ranges that will be encountered in the model domain. This method will be used for 
variables such as the saturation vapor pressures, latent heats, specific heats etc. These 
variables are dependent on temperature alone. Their look-up tables will be generated over 
the temperature range -100°C to 50°C at 0.1 degree intervals. The valué of the variable at 
any given temperature can then be computed by simple linear interpolation between valúes 
picked from the table. 


The following steps form the computational sequence that will be discussed in detail in the following 
sections: 
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1) Evalúate local temperature from potential temperature 

2) Evalúate microphysical parameters 

3) Evalúate local super-saturation 

4) Calcúlate the ice nuclei concentration 

5) Calcúlate the particle size distribution parameters for rain and snow 

6) Calcúlate the production rates of non-precipitating hydrometeors 

7) Adjust mixing ratios of non-precipitating hydrometeors 

8) Calcúlate the production rates of precipitating hydrometeors 

9) Adjust mixing ratios of non-precipitating hydrometeors 

10) Adjust temperature to inelude the release of latent heat 

11) Compute the precipitation (rain and snow) fluxes due to their fall velocities. 


4.3.2 Evalúate Local Temperature from Potential Temperature. 

In OMEGA, the prognostic thermodynamic variable is p0, where p is the density and 9is the potential 
temperature. (p9 in effect is the energy density.) The potential temperature of a parcel of air is defined 
as the temperature that air parcel would have if it is moved adiabatically to a reference pressure level (in 
meteorology the reference pressure level is assumed to be 10 5 Pa, = 1000 mb). The majority of the 
microphysical processes, however, are driven by temperature and not potential temperature. Henee the 
temperature has to be calculated from the potential temperature. The temperature and the potential 
temperature are related via the Poisson equation, 



( _ 

K 

( _ \ 

T = 9 

P 

II 

P 




P,eS J 


where 


(103) 


P Kf =1x10 5 Pa 


and K = — = 0.286 
c„ 


(104) 


4.3.3 Evalúate Microphysical Parameters. 

Several physical parameters are involved in the microphysical calculations. All the microphysical 
processes involving phase change of any type need the respective latent heats. The phase change 
processes which involve the vapor phase also depend upon the local saturation condition. Henee the 
corresponding saturation vapor pressures have to evaluated. Other processes involving heat transfer 
need the specific heats of water, ice, vapor, and air, and parameters such as the diffusivity of vapor are 
required for computing the rate of sublimation or condensation of the precipitating hydrometeors. Once 
the local temperature is defined, these microphysical parameters can be evaluated. 


4.3.3.1 Latent Heats. Latent heats of vaporization and fusión are assumed to vary only with 
temperature, T such that (These expressions are from Pruppacher and Klett (1978).) 


U(T) 


2.50078X10 6 


í T oY 



J/Kg 


(105) 
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where y = 0. 167 + 3.67 X10" 1 7 and T 0 = 273.16 K (= 0°C) and 



3.3369 xlO 5 + 2030.6(7-7 0 ) 
-10.467(7’ -7 0 ) 2 
3.3369x10 5 J/Kg 


J/Kg 


if 7 < 7 0 
if 7 > 7 0 


The latent heat of sublimation is given as 

L s {T)=L v {T)+L f {T) 


( 106 ) 


(107) 


The latent heat of sublimation varíes only very slightly with temperature and can be assumed to be a 
constant at 2.83237xl0 7 J/kg. Henee the following scheme is generally used for calculation of latent 
heats. For temperatures above melting, the latent heat of vaporization is computed using (4.48) and the 
latent heat of fusión is calculated as 


L f =L s -L v (108) 

If temperature is less than melting then the latent heat of fusión is calculated as a function of temperature 
using (4.49) and then the latent heat of vaporization is calculated as 

L v =L 5 -L, (109) 


43.3.2 Specifíc Heats and Other Parameters . The specific heat of water and ice are given by 
Pruppacher and Klett (1978) (in units of J/Kg-K as 


c.(t)- 


5400.0 

4217.8 + 0.3471(7 -7 0 ) 2 
4178.0 + 0.01298(7 - 308.16) 2 
+ 1.591xl0" s (7 -308.16) 4 


if 7 <223.16 
if 223.16 <7 <T n 


2106.1 +7.327(7-7 0 ) 
2106.1 


and C,(t) = 

The dynamic viscosity of air is given by List (1971) as 

¿i d {t)= 1.8325x1 O -5 


if 7 0 < 7 


if 7 < 7 0 

if 7 > 7 n 


416.6 ( 7 Í 5 


7 + 120.0 


296.16 


= 1.49629 x1o -6 

Then the molecular viscosity is simply 

V,(P,T)- 


( T*-5 A 


7 + 120.0 


V 


J 




The thermal conductivity of air is (Wisner et al., 1972) 


( 110 ) 


( 111 ) 


( 112 ) 


( 113 ) 
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Mr)=i4i4.0A>(r) 


-3 


( T 1 - 5 ^ 


= 2.11575x10 . 

T + 120.0 

The diffusivity of water vapor in air is (Pruppacher and Klett, 1978) 

/101325.0 Y T Y 94 


(114) 


D v (P,r) = 2.11x10 


-5 


= 4.0119x10" 


P 

y 1.94 


\ T 0J 


(115) 


The Schmidt number is defined as the ratio of the molecular viscosity to the diffusivity of water vapor in 
air or 


A r St (?.7') = 


y.(p.r) 


-0.56 


D,(P,T) 


= 3.72963 xl0“ 2 /L 


T + 120.0 


(116) 


4.3.4 Evalúate Local Supersaturation. 

The processes of condensation/evaporation and deposition/sublimation are very dependent on the 
amount of water vapor that exists in excess of the saturation limit at that point in space. Henee the local 
saturation vapor pressures (over water and ice surfaces) are calculated. The saturation vapor pressures 
are functions of temperature only. Their variation with temperature is given by the Clausius-Clapeyron 
equation (Hess, 1959). This is somewhat cumbersome to use in a program as it involves other 
temperature dependent quantities such as the latent heats. We will adopt the method used in TASS 
(Proctor, 1987) in which the saturation vapor pressures (over water and ice surfaces) are computed using 
the following polynomial fits to experimental data (Wexler, 1976, 1977): 

e jw (r) = exp(-2991.27297’ _2 -6017.01281"' +17.87643854 

- 0.028354721 T + 0.17838301 x 10~ 4 T 2 

- 0.84150417 x 10 -9 T 3 + 0.44412543 xlO -12 T’ 4 
+ 2.858487 lnT) 

and 

e si {T) = exp(- 5865.36967’"' + 22.241033 + 0.013749042 T 

- 0.34031775 x 10 -4 T 2 + 0.26967687 x 10" 7 T s 
+ 0.691865 llnT) 


The saturation mixing ratios of water vapor over a water surface Q vs , and over an ice surface Q s ¡, can 
then be calculated from e sw and e s ¿, respectively as 



0 _ e swÍQv +£ ) 

*ívs p 

(119) 

and 

Q _ e ÁQv +£ ) 

>¿ 5i p 

(120) 


where e -R/R D - 0.622. 


(117) 


(118) 
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4.3.5 Calcúlate the Ice Nuclei Concentration. 


Ice nuclei (IN) are aerosols on which water vapor can deposit to form ice crystals. Only a fraction of all 
the aerosols present in the atmosphere can serve as ice nuclei. This fraction is determined by the surface 
and crystalline properties of the material. Among all the ice nuclei present only some are active at any 
given temperature. Observations of the number concentration of active ice nuclei as a function of the 
local temperature (Fletcher, 1962) yield 


n 


IN 


10 9 if T < 230.95 

P IN exp{0.6(r M -7)} if 230.95 <T <T, 
0 


M 


if t m <t 


( 121 ) 


where fi IN =0.01 m 3 . The IN concentrations can vary by as much as a factor of 10 at a location and 
temperature (Pruppacher and Klett, 1978). 


4.3.6 Calcúlate A’s and Fall Velocities for Rain, and Snow. 

The intercepts of the inverse exponential distributions describing rain, and snow are defined as 


and 


N 0R = 2.25 xlO 7 


f Q l l ^ 0375 

kP J 


( 122 ) 


N os = 5.5 x 10 6 exp[- 0.088(T - T M )] 

The slopes of the inverse exponential distributions can now be calculated as 




pQr 


\0.25 


and A s = 


pQs 




25 


\^7iN os S s 


(123) 


7lN 0R S w 

< 

where 8 W and 8 S are the densities of water, and snow respectively given as. 

¿ w = 1000. Kg/m 3 and S s =100. Kg/m 3 

The precipitating hydrometeors have terminal fall velocities, which are calculated from a single particle 
fall velocity averaged over the respective distribution in a mass-weighted manner. The mass-weighted 
fall velocities for rain, and snow are calculated as follows: 


= (- 0.267 + 2.06 x 10 4 A R - 2.045 x 10 7 A 2 * + 9.06 x 10 9 A 3 * í- 


L2 

P 


x0.4 


(124) 


and 




1.1 


/ \0.5 

' 1.2 ' 


if T < T, 


2.5 


KP J 
'12** 


M 


(125) 


J 


if T > T t 


M 


4.3.7 Production Rates of Non-Precipitating Hydrometeors. 

With the local thermodynamic State and key microphysical parameters computed, it is now possible to 
evalúate the production/loss of the hydrometeors. Because their smaller size increases the surface area 
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to volume ratio, the non-precipitating hydrometeors respond most quickly to changes in their 
surroundings. For this reason we evalúate their changes first. In the following discussion of production 
terms, we will use the following terminology. Q will denote the mixing ratios, AQ the change in Q, and 
P the rate of change of mixing ratios with time. The subscripts of AQ and P will indícate the 
microphysical process in volved. The non-precipitating hydrometeors are comprised of the cloud 
droplets and the ice crystals. Even though they have finite fall velocities, they are small enough to be 
ignored, and for all practical purposes these partióles can be assumed to follow the airstream. For the 
sake of computational efficiency these calculations are performed in a model layer only if at least one of 
the participant species is significantly non-zero in some cell in that layer. 

4.3.7.1 Freezins of Cloud Droplets. Freezing of water is initiated by ice nuclei that are immersed in 
water, or come into contact with the water surface. Large water drops generally contain enough 
impurities that they freeze at or near 0°C (Tm). Small droplets, however, can exist in liquid phase at 
temperatures far below 0°C, if they do not contain or come into contact with any ice nuclei. Even small 
droplets however will spontaneously freeze at about —40°C. In nature, supercooled liquid is usually 
found in very vigorous thunderstorms, where the droplets are formed in a rapidly cooling air parcel in 
the strong updrafts. In this versión of OMEGA, however, no supercooled liquid is allowed to exist. 
(Supercooled liquid and melting solids will be allowed to exist in a future versión.) Therefore, if 
T < Tjy[, all cloud droplets in that computational cell are frozen to produce ice crystals as shown in the 
following relation: 


Q¡ — Q¡ + Qc 

A Qcfzn ~ Qc 

Q c = o. 


if T <7) 


M 


(126) 


4.3.7.2 Meltins of Ice Crystals. If T > T M all ice crystals in that cell are melted to form cloud droplets 
by the following: 


Qc ~ Qc Q¡ 
A Qmn = Q, 

Q¡ = 0 . 


if T>T h 


M 


(127) 


4.3.7.3 If T> TCompute Condemation/Evavoration . The amount of cloud water vapor condensed 
into liquid depends on the water vapor excess (over saturation) that is present at a grid point. The excess 
vapor is assumed to condense into cloud droplets (or undergo deposition to form ice crystals.) This 
process is assumed to occur instantaneously; thus supersaturated conditions cannot exist in the model 
domain. (This process will be replaced by a condensation rate equation later, which will then allow for 
supersaturated conditions to exist.) 


The condensation adjustment equation can be derived by considering a closed volume of moist air. If 
supersaturation exists, the excess water vapor will start to condense. As condensation proceeds, latent 
heat of vaporization is released, increasing the temperature and henee the saturation vapor pressure of 
air. Thus the initial water vapor excess AQ V is used not only for condensation, but also for the increase 
in vapor capacity of the air. Let AQ¡, and AQ 2 denote the amount of water vapor condensed, and the 
amount of extra vapor the air can hold due to the increase in temperature such that the total AQ is given 
by 
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AQ = AQ ] +AQ 2 


( 128 ) 


The increase in temperature A T due to the release in latent heat when AQ¡ vapor is condensed is 

(i 29) 


A T = 1 

c„ 


This increase in temperature increases the saturation vapor pressure as given by the Clausius-Clapeyron 
equation, 


de^ = L±dT 
e, R T 2 


(130) 


In an isobaric process (we assume that in each grid cell the pressure does not change during the 
microphysical adjustments) such that 

dQ vs _ de s 


Q v¡ 


(131) 


henee, 


AQ- 


L V AT L v 2 Q sv AQ l 

n rr.? n m 9 


R V T 2 c p R v T 2 


(132) 


Define a fraction r such that AQi = r A Q. Then 

A<2, 


i-i 


r = M- 


A Q AQ ] + AQ 2 ¡ [ 


=Afi Ae 1+ ^^ - 
C P RJ 


= 1 + 


WQ„ 

c p R v T 2 


(133) 


Thus the amount of water vapor to be condensed (or the máximum amount of cloud water that can be 
evaporated) can be written as 

a-Ov 


A Qccnd ~~ ^A Q — 


i , l?q„ 

c p R v T 2 


(134) 


4.3.7.4 Initiation of Ice Crystals . Ice crystals can be produced as vapor deposits directly onto ice nuclei 
(deposition nucleation). The production rate of ice crystals via this process thus depends on the ice 
nuclei number concentration, which in tum depends on the local temperature. 

In OMEGA cloud ice is initiated at temperatures less than -16°C; at higher temperatures, the rate is near 
zero. The initiated crystals are assumed to ha ve a máximum diameter of 12.9 pm and a máximum mass 
of 10~ 12 kg. The amount of ice initiated by this process is limited by the assumed máximum ice crystal 
mass (and the IN concentration) and by the available supersaturation (with respect to an ice surface). 

The máximum ice mixing ratio attained by growing all crystals to 12.9 ¡Jm is 

AQ n = 10 12 ” w (135) 

P 
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The máximum ice mixing ratio that can be generated by depositional growth driven by the 
supersaturation can be derived in a manner analogous to the condensation adjustment (Equation 97) by 
replacing L v by L s , and Q vs by Q s ¿. Thus, 


A O = -^ v _ Qá— 

12 L 2 Q. 
1+ 


c p R v T 


(136) 


and 


^Qicwv i — min(A<2/i , AQ !2 ) 


(137) 


4.3.7.5 If T <Tm Compute Deposition/Sublimation . Existing ice crystals grow by deposition if, after 
ice crystal initiation, excess vapor (over ice saturation) still exists. The depositional growth rate of ice 
crystals is given by 

dgj _ 65 2 ¡ QiP 8»~8.| 

3t ' 1¡n„ Q.A’+B' 

where 

« L ( L 
k T T{R v T J 

and fi*=— (140) 

PD V 

where D v is the diffusivity of water vapor given by Equation 115. 

This rate has to be limited by the máximum cloud ice that can grow by deposition under the current 
saturation conditions. This is given by A Qj 2 ', henee, 

Pidep - m ¡ n í-jp- (MI) 

at At J 

If {Qv-Qsi) is negative then sublimation occurs. (The sublimation model will need refinement to better 
define cirrus clouds. Observations suggest that ice crystals from cirrus can fall large distances through 
sub-saturated air.) 

4.3.8 Adjust Mixing Ratios of Non-Precipitating Hydrometeors. 

Before proceeding to compute the production/loss of precipitating hydrometeors, the mixing ratios of 
cloud droplets, ice crystals and water vapor are adjusted. 

Qc ~ Qc A Qccnd 

Q,=Q,+ R¡dep A^ A2/cwi (142) 

Q v —Q v ~ A Qccnd ~ ?idep A* ~ ^Q/cwv i 
The temperature adjustment is determined from 

c pAT = (AQcfzn ~ AQimlt )A/ + AQ CCND L V + (A<2 /CW1 + P IDEP At)L s (143) 


(138) 


(139) 
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4.3.9 Production Rates of Precipitating Hydrometeors. 

Having completed the adjustments of the non-precipitating hydrometeors, we are now able to evalúate 
the production/loss of the larger precipitating hydrometeors. 


4.3.9.1 Production o f Rain bx Meltins ofSnow. In the current versión of OMEGA, all snow melts 
instantaneously to rain above 0°C (Tm) such that 


Qr — Qr + Qs 
A Qsmlt ~ Qs 
Qs =o. 


tf t>t ; 


M 


(144) 


4.3.9.2 Production ofSnow bv the Freezine of Raindrovs . This reverse process to the production of 
rain by melting of snow assumes that all rain freezes instantaneously to snow below 0°C where 


Qs -Qs + Qr 

^■Qrfzn = Qr 
Qr =0. 


if T <T t 


M 
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4.3.9.3 Autoconversion of Cloud Droplets to Form Rain. Even though we have assumed the cloud 
droplets to be monodisperse, in reality they exist in all sizes (small of course). They undergo 
interparticle collisions and coalesce to form larger droplets or rain drops. This stochastic process of self 
collection and growth to form rain is termed autoconversion. The rate of production of rain from cloud 
droplets via autoconversion is given by (Berry and Reinhardt, 1974a and 1974b) 


Pracv —7.26x10 


10 


20 


0.38 


r e ' -0.4 


pQ c 


í o \X 


0.38 


r c 10 6 -7.5 
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where r, = 


' l£üc_ ^ 

4 XNrS, 


V 


c v w 


is the mean cloud droplet radius, a is the dispersión coefficient of the cloud 


droplet distribution, and, Nc the cloud condensation nuclei (CCN) concentration. 


4.3.9.4 Collection o f Cloud Droplets b\ Raindrov . Raindrops collect cloud droplets as the raindrops 
have a much higher fall speed than the cloud droplets. The collection efficiency is a function of the size 
of cloud droplets. Large droplets will have a high collection efficiency, whereas the very small droplets 
will flow around the raindrop and thus have a collection efficiency of zero. The production rate due to 
this collection mechanism is 


Frac ~ 


15;r 

~38"" 


^■r N 0 rW r Q c E RC d 


( 147 ) 


where E RCD is computed from a polynomial fit to experimental data (Proctor, 1987) given by 
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o 


if r c <\2fim 


- 0.27544 + 0.26249 x 10 6 r c 

- 1.8896 xlO'V +4.4626 xlO'V 

1 


if í.2jum <r c < 20ftm 
if 20 fjm < r c 


(148) 


4.3-9.5 Condensation Evaporation ofRaindrovs. Raindrops evapórate in sub-saturated conditions, if 
enough cloud droplets are not present to alleviate the saturation déficit. (The smaller cloud droplets are 
allowed to evapórate first as the saturation vapor pressure over a curved surface with smaller radius of 
curvature is higher than over a surface with larger radius of curvature.) Also if the environment is 
supersaturated, vapor can condense onto the rain drops. The following expression takes into account the 
ventilation effect of the falling raindrop 



4.3.9.6 Autoconversion oflce Crvstals to Snow. Cloud ice conversión to snow occurs when the 
average crystal mass exceeds 9.4 xlO -10 Kg. This assumes that there is a stochastic máximum ice 
crystal size. The excess ice crystal mass is assumed to contribute to the snow field. This process is 
analogous to the production of rain by autoconversion of cloud droplets. 

Q\ ma x (th e máximum cloud ice concentration) = 9.4 x 10 -10 n, N (150) 

íw=(a-a, j-J- asi) 

A t 

4.3.9.7 Collection oflce Crvstals bv Snow. Just as raindrops collect cloud droplets by virtue of their 

difference in slip velocity, snowflakes collect ice crystals. The production rate of snow due to this 
collection of ice crystals is given by 

P SACI =Q.5nN os W s A s 3 Q'E sa (152) 

where E sa , the collection of efficiency of snow for ice crystals, is a function of temperature such that 

£ sc ,=exp(0.38(r-r„)) (153) 

4.3.9.8 Deposition/Sublimation of Snow . The deposition onto (or sublimation of) snow is calculated in 
a manner similar to condensation on (or evaporation of) raindrops. The following formulation also takes 
into account the ventilation factor of falling snow: 
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( 154 ) 


4.3.10 Adjust Mixing Ratios of Hydrometeors. 

With all the production/loss terms evaluated, the mixing ratios of the hydrometeors and water vapor are 
adjusted by 

Qr = Qr **■ (Pracv ■*" Prcnd "*■ Frac W 

Qs ~ Qs + i^SACV "*■ ^SDEP ■*" ^SACI W 

Qc ~ Qc ~ ( Pracv "b Frac (155) 
Ql=Ql~ ( PsACV + PsaCI W 
Q v = Q, ~ {.PrcND ■b PsDEP W 

4.3.11 Adjust Temperature to Inelude Latent Heat Release. 

The final step is to adjust the temperature in the cell to account for the latent heat released or absorbed 
during condensation/evaporation, deposition/sublimation, and freezing/melting involving the 
precipitating hydrometeors. The total latent heat released is balanced by the change in sensible heat in 
the following formulation: 

CpAT = (áQ rfzn — AQ smlt )L f + P RCND AtL v + P SDFP AtL s (156) 

This equation determines the temperature difference due to latent heat release. After the temperature is 
updated, the new valúes of potential temperature are calculated using the Poisson equation. 

4.3.12 Compute Precipitation Fluxes. 

The fluxes due to precipitation are straight up and down, through the top and bottom faces of each cell. 
Even though the precipitation contribution can be accounted for by adding the fall velocities to the 
advection velocities, it is computationally more efficient to calcúlate it separately. The mixing ratio 
changes to rain and snow can be calculated as follows. 

and % = i (,57) 

at p az atpáz 

where z denotes the local vertical direction (not to be confused with z the vertical coordínate of the 
Cartesian model domain). 

4.3.13 Convective Parameterization. 

Deep cumulus convection is associated with a deep conditionally unstable layer and the presence of 
large-scale convergence. The first condition makes it possible for huge cumuli to penétrate into the 
upper troposphere and the Iower stratosphere, while the second condition provides a general lifting 
mechanism to trigger the convective instability. Therefore, deep cumulus convection always takes place 
in regions of conditionally unstable stratification and low-level convergence of moisture, and in such 
regions this process causes the release of latent heat. 
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The release of latent heat by cumulus convection is one of the most important short temí effects on the 
atmospheric dynamics. In addition, the water droplets within the clouds play a prominent role in the 
reflection, absorption, and emission of solar and terrestrial radiation. However, due to the large 
horizontal space scale and time steps of mesoscale models the vertical transpon of heat, moisture and 
momentum by deep cumulus convection may not be adequately resolved. While the goal of OMEGA is 
to try to explicitly resol ve large areas of convection, there will always be regions which are not 
sufficiently resolved. To circumvent this problem a versión of cumulus parameterization that was 
originally proposed by Kuo (1965,1974) and later modified by Anthes (1977a) is incorporated to 
account for the effect of subgrid scale deep cumulus convection on the local environment. This 
convective parameterization only applies to those regions which are convectively unstable to deep 
penetrative convection and in which the total horizontal moisture convergence exceeds a critical valué. 
As the cumulus parameterization is a mechanism to account for subgrid convection in large cells, it is 
usually applied only to cells that large in area (0(100 km 2 )). In OMEGA we are exploring a concept by 
which we can smoothly transition from regions where no convective parameterization is applied to 
regions where it is applied. This will be achieved by including a factor in the cumulus adjustment 
scheme which can vary between 0 and 1 depending on cell area. However, the explicit cloud 
microphysics is active over the whole domain. 


The total rate of moisture accession, M t , per unit horizontal area is given by the convergence of moisture 
in the column of atmosphere above the unit area plus the surface evaporation, thus 


M 


= - / P (V„ yq)dz+S e ¡p 0 C o V ( g¡ - q a )] 


(158) 


where zt is the model top, q the cell averaged mixing ratio, q s the saturation valué corresponding to sea- 
surface temperature, q a the air valué just over the sea, C D the drag coefficient, p the density, p 0 air 
density in the surface layer, V 0 wind speed at the surface layer, and S c is an indicator such that 


S = 


c 



overland 
over water 


(159) 


Henee, evaporation within the surface layer is taken into account only over the sea when w 0 > 1 cms' 1 
96 

and —— < 0 , where w 0 is the vertical velocity in the surface layer, and 0 e is the equivalent potential 
oz 


temperature. It is assumed that the direct evaporation effects do not contribute to the generation of deep 
cumulus convection over land. The drag coefficient is determined using the surface layer similarity 
relationships as follows: 
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where /? is a constant and set to 1.35, k von Karman’s constant, and X ¥ H the integral functions (see 
PBL parameterization), zq the surface roughness height. 
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The moisture M t is then used to make cloud columns (T c , q c ) from environmental air (T ,q). Part of 
the moisture (Ai,) will condense and precipítate, raising the temperature from T to T c and the remainder 
(M 2 ) will increase the mixing ratio of the cloud column from q to q c , that is 


c 


M = -*- 
1 L 


jp(T-T)dz 


c Z 


cb 


(161) 


and 



(162) 


where Cp is the specific heat at constant pressure, L c latent heat of condensation, z ch and z c , the bottom 
and top of the cloud, respectively. Henee the total amount of water vapor needed to create the cloud 
over a unit area is then given by the sum of the two integráis 


M 


= M, +M. 


ct 

=1 


cb 



dz. 


(163) 


The rate of cloud production per unit time, C (or fractional area), is assumed to be proportional to the 
convergence of water vapor plus evaporation for the column divided by the amount of water vapor 
necessary to produce the model cloud as follows: 

C = TT- (164) 

M c 

Then, the product C At gives the total production of cloud air during time At. 


In 1974, Kuo also introduced a división of the convergence of moisture into a fraction bM, of the total 
which increases the humidity of the air, and a fraction (l-b)M¡, which is condensed and precipitates as 
rain or is carried away with the latent heat warming the air (Kuo, 1974). A proper evaluation of the 
parameter b is important. Anthes (1977a) proposed that b depends on the mean relative humidity in the 
air column in such a way that the moistening associated with cumulus convection is strong when the 
atmosphere is dry and the moistening is weak when the atmosphere is wet so that 


(i-rao 

1 


if RH > RH C 
if RH < RH C 


(165) 


where RH is the mean environmental relative humidity in the air column, RH C is the critical relative 
humidity, and n a positive exponent. Henee, when RH is less than RH C , only environmental moistening 
takes place. Anthes (1977b) chose valúes of n=l and RH C = 0.5 for the axisymmetric hurricane model. 
However, a number of studies have indicated that the b valué computed with n=l and RH C = 0.5 is too 
large and consequently convective heating, drying, and rainfall rates are underestimated. In a study on 
semi-prognostic tests of Kuo-type schemes in an extratropical convective system, Kuo and Anthes 
(1984) found that the simulated rainfall rate agrees best with observations when n is between 2 and 3, 
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and RH C is between 0.25 and 0.50. During the Asían summer monsoon periods, Someshwar Das et al. 
(1988) found that the best results are obtained for valúes of n between 3 and 5, and RH C between 0 and 
0.25. In the ECMWF numerical forecasting model, valúes of n=3 and RH C = 0 are used (Tiedtke 1986). 
As a compromise between these valúes, in the OMEGA model n and RH C are specified as 3 and 0.5, 
respectively. 

Now assume that the large-scale forecast is made for a time step A t giving a temperature T e without the 
effects of cumulus convection. Then, applying the parameterization scheme gives a temperature 
T(t + At) suchthat 

r(í+Aí)=7;+(i-¿)cAí(r-rj, ü66) 


and the corresponding equation for mixing ratio q(t + Ai) after application of parameterization is 

q(t + At) = q e +bC At(q -q e ). 

' C ' 


(167) 


The precipitation from the subgrid cloud scale is that parí of the moisture used in warming the air from 
T e to T. Therefore, the rate of precipitation per unit time per unit volume in the column is given by 


pre q-t)cc,(r-r.) 

el J, J z 1 


dz dt. 


(168) 


* On evaluating the subgrid scale cloud effects, cumulus convective processes are activated when the 
vertically integrated moisture convergence rate exceeds a critical valué of 10~ 5 kgm ’ V 1 (Anthes 1977b). 
The cloud base z cb and cloud top z a need to be determined. The cloud base is determined by the lifting 
condensation level (LCL) of each air parcel. The temperature T c and mixing ratio q c inside the cloud are 
considered as those on a moist adiabatic passing through the lifting condensation level. Note that q c is 
the saturation mixing ratio at the temperature T c . The cloud top is specified as the height where cloud 
temperature T c becomes less than the environment temperature T e . 


4.4 Radiation. 


The radiative source-sink term in the conservation of energy relation can be written as 



(169) 


where the terms on the right hand side represent the temperature change resulting from longwave and 
shortwave radiative divergence flux in the vertical direction. The divergence of radiative energy in the 
horizontal direction is neglected, since its variation is much larger in the vertical direction on the 
mesoscale. The methods of parameterizing this vertical flux takes into account the absorption of 
shortwave radiation by water vapor and the longwave energy emitted by water vapor and carbón dioxide. 
It is essentially similar to the one used by (Mahrer and Pielke,1977). Because of the separation of 
wavelength in the atmospheric radiation spectra it is convenient to develop sepárate parameterizations 
for long and short wavelengths. 
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4.4.1 Shortwave Radiative Flux. 


The diuraal variation of the solar flux on a horizontal surface at the top of the atmosphere is computed 


from 




S = S 0 cosZ, 

(170) 

with 




eos Z = eos \¡/ eos 8 eos H + sin y sin 8 

(171) 


where S 0 is the solar constant, Z the zenith angle, y/the latitude, ¿the solar declination, a function of 
Julián day, and H the solar hour angle. Assuming that shortwave absorption in the atmosphere is only 
due to water vapor, the heating of the atmosphere by this radiation is then given by 




Í0.231 S « CMZ 
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where r(z) is the optical path length of water vapor above the layer z and is given by 

r(z) = \'° p pqdz . (173) 

This is the shortwave flux formulation used in the OMEGA model. 


4.4.2 Longwave Radiative Flux. 

The parameterization of the long wave radiative flux in atmospheric models is typically treated as a 
function of the normal optical thickness which when integrated over all wavelengths is represented by 
the broadband emissivity, e . In clear or cloudy air this emissivity is dominated by the water content of 
the air. Water vapor and carbón dioxide (CO 2 ) are considered as emitters of longwave radiation. The 
path length for water vapor (Ar ; ) expressed in units of g cm' 2 is computed for each vertical layer by 

ÍP -P) 

A 0 =Pj<lj(Zjn-Zj) = -—- -qj ; (174) 

o 

and the path length for carbón dioxide (A Cj) expressed in units of millibars is 

A cj = -0.4148239(p /+1 - Pj ). (175) 

After these increments are obtained, they are summed up from the fírst level to the ith level to give the 
total path length as follows 

J/+i = ¿Ar- and c I+1 = ¿Ac- (176) 

j= 1 7=1 

where ri = ci =0 at the surface. 


The emissivity for water vapor was derived from data in Kuhn (1963) and is given in Jacobs et al. 
(1974)by 
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£ r ( i J)= 


0.113 log 10 (1+12.63F) log I0 < -4 

0.104log )0 r+ 0.440 -4<log 10 r<-3 

0.121 log 10 r +0.491 -3<log ]0 r <-1.5 

0.146 log 10 r + 0.527 —1.5 < log )0 r < -1.0 

O.1611og 10 r + 0.542 -1.0<log 10 F<0 

0.136 log 10 r + 0.542 log 10 r>0 


(177) 


where r = j/- — r ; j is the optical path length between the ith and jth level. The emissivity function for 
carbón dioxide is given by Kondratyev (1969) as 

£ CÜ2 (i, j) = 0.1851 -exp(-0.3919c ¿ -° 4 ]. (178) 

Total emissivity for each depth between level i and level j is then given by 

£(U j ) = £ r (U j ) + £ C o 2 & j ) • ( 179 ) 

The downward and upward radiative fluxes at a level N can be computed using the above emissivity 
functions as 


tOP~ 1 . , 

R i (N) = X-fe, + t; \e{N, j + 1)- £ {N, j)]+ ¿r; op [1 - e{N,top)] 


(180) 


and R t (N) = T , 4 +1 + T¡ )[e{N, j) - e{N, j + 1)] + < [l - £ {N ,\)]. 

where Te and T top are the temperatures at the ground level and model top, respectively and a is the 
Stefan-Boltzman constant. Thus the radiative cooling at each level N except the ground level is 
computed as 

1 [Rl (N + V-Rl (N) + Rl(N)-Rl(N + \)] 
j LW pc p z(N + \)-z(N ) 

where LW denotes the longwave radiation and z is the height. 


(181) 


(182) 


Without simplification the computation of radiative transfer is computationally expensive. Sasamori’s 
technique (1972) which assumes an isothermal atmosphere for radiative transfer that simplifies the 
computing procedure. After this simplification, the temperature change resulting from longwave 
radiative flux divergence at each level N is computed as 


(c¡T¿ - oT¿\ £ (N + 1,1) - £ (N ,1)]+ {oT* - oT„ \s(N + l,top) - £ (N,top)] 


z(N + l)-z(N) 


(183) 


This is the method used in OMEGA. 
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Section 5 


Grid Generation and Adaptation 


5.1 Grid Definition. 

As mentioned in Section 2, the OMEGA model uses a structured-unstructured grid consisting of vertical 
columns of prisms: five-sided polyhedral having triangular tops and bottoms and quadrilateral sides. 

The grid is unstructured in the horizontal direction (there is no imposed ordering of triangles) and 
structured in the vertical (indexing in the vertical direction is sequential and is predetermined). The top 
and bottom faces of the prisms have the same radially projected area, but each triangle is allowed to tilt 
independently. The surface layer, which is generated first, conforms to the terrain. This layer is then 
projected radially outward from the earth’s surface, resulting in layers of triangular grids that form the 
tops and bottoms of vertical columns of prisms. The terrain roughness is smoothed over a specified 
number of layers, resulting in a fíat (i.e. horizontal) top layer. The complete grid can be described by 
describing the triangles forming the surface layer and the altitudes of vértices in the upper layers. 

The description of the unstructured grid requires pairs of arrays containing information about each 
vertex, edge (line segment connecting two vértices), and triangle. Each pair consists of a real array 
containing spatial information (vertex coordinates, edge lengths, and triangle areas), and an integer array 
containing connectivity and boundary information. Because the grid is structured in the vertical 
direction, spatial and connectivity information need only be stored for the surface layer. Within the 
OMEGA code, adjustments to lengths and areas are made based on radial distance from the earth’s 
surface. 

The primary benefit of the unstructured grid over a conventional structured grid lies in the ability of 
smooth transition from high resolution where needed to low resolution elsewhere. Currently, the 
OMEGA grid generator automatically adapts a grid to topological features that affect weather, such as 
teiTain roughness and land/water boundaries, producing smaller triangles in these areas. Over oceans, 
where variations in the horizontal are small, larger triangles are created. By substituting another dataset, 
the grid can be adapted to any geographical property. For example, by interfacing with a file containing 
the current meteorological conditions, the grid can be adapted to any field quantity such as pressure or 
temperature, producing higher resolution in areas of large gradients, thus allowing higher resolution in 
regions of frontal activity or convection. This latter capability has been developed and tested, but is not 
interfaced with the current production versión of the OMEGA modeling system. 

An example of the flexibility of the OMEGA grid is shown in Figure 17. This figure shows a grid 
generated for México in which the grid was adapted to the underlying topography, the land/water 
boundary, and to the initial weather conditions. The synoptic situation chosen was Typhoon Linda. In 
this example, we have broken the grid generation process into different steps for iílustration. We have 
shown a grid generated by adapting to gradients in elevation (refining the grid in mountainous areas), 
gradients in the land/water Índex (refining the grid in Coastal areas), adapting to typhoon Linda (refining 
the grid to the initial weather), and, finally, to all of these criteria. 
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(a) Adapting to topography. (b) Adapting to land/water boundary. 



(c) Adapting to weather (Typhoon Linda). (d) The resulting grid. 


Figure 17. Using an appropriately chosen set of criteria, OMEGA can generate a grid that 
captures the important physical features. 

In addition to the automatic adaptation to geographical features, the grid can be adapted to any user- 
defined geographical area, such as a theater of operation, by the creation of up to ninety-nine rectangular 
subdomains on which higher resolution can be specified. Within each subdomain, the adaptation is 
govemed by the user-specified minimum and máximum resolutions. If these limits differ by a factor of 
two or less, the resulting grid will consist of triangles that appear uniform to the eye. Subdomains can 
overlap to create a high-resolution región of almost any shape. They can also be layered to produce a 
high-resolution area within an intermediate-resolution región. In addition to subdomains, the user can 
also define a point on the grid and a radius of influence around it. The grid will automatically refine 
within the región of influence. A smoothing process is performed on the resolution specifications to 
avoid high expansión ratios at the boundaries of subdomains. 
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5.2 Grid Generation. 


The OMEGA grid is constructed using several routines that add, subtract, and relocate vértices to satisfy 
preset specifications. The grid generation methodology consists of the following steps: 

1. Create files in the working directory containing the terrain and land/water data, culling 
information from the global databases. 

2. Generate an initial coarse, two-dimensional (terrain following) grid of triangles. 

3. Refine the grid through a series of iterations, calculating vertex elevations and cell properties 
(currently land or water) by interpolation using the database after each iteration. 

4. Create the vertical levels, smoothing the roughness gradually in layers above the surface of the 
earth. 

5. Generate the surface characteristics file. 

These steps will be discussed in the following paragraphs. 

5.2.1 Data Files. 

Two sets of global terrain and land/water data are currently available for the creation of grids: a fine 
resolution database and a coarse resolution database. The coarse database currently in use consists of 
two files. The terrain file contains data for the entire globe at a resolution of five minutes, or about 9.26 
km at the equator, while the land/water file contains data at a resolution of ten minutes. The fine 
resolution databases of both terrain and land/water, are divided into sets of files each of which contains 
data for a five degrees by five degrees segment of the globe. The resolution of the fine data files is given 
in Table 2. Any combination of high and low-resolution terrain and land/water databases can be used for 
the primary domain and subdomains of a grid. While high-resolution databases will almost always result 
in a grid that captures geographical features more accurately, using a low-resolution database saves on 
disk space and greatly speeds up the grid generation process. 

Table 2: High-resolution OMEGA terrain and land/water databases. 


Latitude (N or S) 

Resolution (NS x EW) 

80-90° 

30 x 180” 

0.93 x 0.97 km 

75-80° 

30x 120” 

0.93 x 0.96 km 

70-75° 

30 x 90” 

0.93 x 0.95 km 

50-70° 

30 x 60” 

0.93 x 1.19 km 

0-50° 

30 x 30” 

0.93 x 0.93 km 


The first task of the grid generator is to lócate and read the data files. Data that is needed for the case is 
temporarily stored one degree of latitude at a time, in work arrays. This data is then written to files 
named grdnnlan.dat and grdnnter.dat, where nn is the domain or subdomain number. These files are 
read repeatedly, and the data used for interpolation, throughout the grid generation process. 

The low-resolution databases are the default databases. If the grid generator fails to lócate a fine 
resolution file for any domain, the database for that domain will revert to the coarse database. Since 
terrain and land/water datasets are created independently, coarse data can be used for terrain definition 
while fine data is used for land/water definition, and vice-versa. 

Because the land/water Índex for each cell takes the valué zero or one, the land/water data is stored as a 
formatted file using an f3.0 format, resulting in a significant savings in disk space. The terrain file is 
written as a binary file, which speeds execution. 
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(a) Before. 


(b) After triangularization. 


Figure 18. Initial rectangular grid. 


5.2.2 Initial Grid. 

The boundaries of the domain are input as longitude and latitude, while the desired resolution is 
specified in kilometers. For construction of the initial rectangular mesh, meters per degree of longitude 
is determined at the center of the domain. Meters per degree of latitude is globally constant. Aside from 
the variation from north to south, the spacing of the vértices in the initial rectangular grid is equal to the 
mínimum resolution requested for the primary domain (Figure 18). Since the coordinates of each vertex 
are expressed in longitude and latitude, the grid is rectangular in latitude/longitude space. 

Starting with a coarse grid and refining minimizes the work done in generating the grid, since the grid 
never contains more cells than are necessary. To prevent the final grid from being too “regular”, 
randomization is produced in the initial grid by randomly selecting a subset of edges and forcing 
reconnection. The reconnection process is described below. 

5.2.3 Refinement of the Grid. 

After generation of the initial grid, the process of refinement begins. The processes used to refine the 
grid were described briefly in Section 2; in this section we discuss more fully the mechanics of the grid 
refinement. The refinement process begins with the calculation of the area of each triangle, which is 
then compared to the areas of equilateral triangles having edge lengths equal to the mínimum and 
máximum resolutions specified in the input file. A triangle is “flagged” for refinement if trisection will 
not result in new cells having areas less than or equal to the area of an equilateral triangle with the 
mínimum allowable edge length, and the cell satisfies one or more of the following: 

• The area of the triangle exceeds the area of an equilateral triangle with edges equal to the 
máximum edge length specified for the domain or subdomain in which the centroid is located. 

• The elevation gradient between the triangle (the average elevation of the three vértices, since 
elevation is a vertex-based quantity) and its neighbors exceeds a preset tolerance, currently set 
to 1.5 percent of the variation in elevation throughout the domain. 

• The cell is on a land-water boundary. The land/water índex is a cell-centered valué, and is set 
to zero or one; the whole cell is defined as either land or water, depending on the location of 
the centroid. If the land/water indices of any of the cell's neighbors (the three cells that share 





edges with the target cell) differ from that of the target cell, the target cell is on a land-water 
boundary. 

• The area of the cell is greater than four times the area of its smallest neighbor. 

Refinement of a triangle is accomplished by adding a vertex at the centroid, if the triangle is in the 
interior of the grid, or by bisecting the boundary edge, if the triangle is on the boundary. A cell is 
flagged for recombination (the vertex opposite from the longest edge will be deleted from the grid) if it 
satisfies any or all of the following: 

• The area of the cell is less than the area of an equilateral triangle with edges equal to the 
minimum edge length specified by the user for the domain or subdomain in which the centroid 
is located. 

• The cell is water surrounded by water, and combining three cells of the same size into one will 
not result in a cell whose area exceeds the máximum allowed. 

• The cell has aspect ratio less than 0.33. The aspect ratio is defined as the height of the triangle 
divided by the length of the longest edge. 

The criteria for both refinement and recombination are “weighted”, that is, certain criteria will forcé one 
or the other action to occur regardless of any other properties of a triangle. These criteria are chosen 
with the OMEGA solver in mind. Triangles that are too large or too small are forced to trisect or 
recombine to prevent adverse effects on the time-step, and triangles with poor aspect ratios are forced to 
recombine to avoid instabilities in the solution. Comer triangles are bisected to avoid triangles having 
two boundary edges. 

Figure 19 shows the grid in Figure 18 after two triangles, one in the interior and one on the boundary, 
have been refined. Adding and deleting vértices results in a grid with unequal edge lengths and new 
triangles that are frequently obtuse. These irregularities in the grid may cause instabilities to occur 
during the simulation; to avoid this the grid is relaxed. Relaxation entails relocating each vertex to a 
point equidistant from its connecting neighbors. This stage is shown in Figure 19. 




(a) Initial rectangular grid 
with two vértices added. 


(b) After triangles are refined the 
vértices are relaxed (open circles). 


Figure 19. Grid refinement. 
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Figure 20. Grid after reconnection. The deleted edges are shown as 
dashed lines and the new edges are shown as thin lines. 

Following relaxation, triangle shape is optimized through the process of reconnection. During 
reconnection, two neighboring triangles are joined to form a quadrilateral. One diagonal of the 
quadrilateral is the shared edge. These triangles are compared with the triangles that would result if the 
other diagonal of the quadrilateral were the shared edge. If the resulting triangles are more nearly 
equilateral, the original shared edge is deleted from the grid, and the new edge added. Figure 20 shows 
the grid in Figure 19 after reconnection. 

An iteration is complete when all of the flagged cells have been altered. Since grid refinement has 
resulted in the creation and relocation of both vértices and cell centroids, interpolation to the underlying 
databases, both terrain and land/water, is repeated. Bilinear interpolation is used to calcúlate both 
elevation and the land/water Índex. A valué between zero and one will result for the land/water Índex, 
which is then rounded to zero or one. 

This process of refinement and interpolation is repeated until the grid is determined to have converged. 
Convergence is declared when two consecutive iterations result in changes to fewer than 95 and 99 
percent of the cells, respectively. The number of iterations required depends largely on the ratio between 
the mínimum and máximum resolution desired, since the initial grid is based on the mínimum 
resolution, and numerous trisections of cells may be required before the máximum resolution is attained. 
Our experience has shown that four iterations are typical when the mínimum and máximum resolutions. 
differ by a factor of two or less, while eight or nine iterations may be required when the resolution ratio 
is five or more. 

A final step prior to the creation of the vertical levels consists of bisecting any comer triangles that have 
two boundary edges. This process has been found to improve the solution at the comers of the domain 
(See Figure 20). 

5.2.4 Vertical Levels. 

The final step in grid generation is the creation of the vertical levels. The OMEGA grid is constructed in 
two parts in order to provide high resolution near the surface, where it is required, without requiring that 
same resolution at altitude. The lower levels are considered “stretchable”, while the upper levels 
represent spherical surfaces. The user is required to specify both the total number of levels, N, as well as 
the number of stretchable levels, N s . 
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Figure 21. The vertical levels of the OMEGA grid are constructed by 

extending radiáis through each of the surface vértices. A set of 
stretchable layers exists near the surface with constant thickness 
layers above. 

In normal operation, the user specifies the height of the first layer of triangles, Az 0 , above the surface of 
the earth along with the máximum geometric factor, a, nxa , by which the heights of the prisms will 
increase in the radial direction in the stretchable layers. The layers are constructed by taking each 
surface vertex and constructing an array of vértices along the radial through that vertex. The resulting 
vértices are then connected to form the OMEGA grid volumes. 

The height of each of the vértices in the interior of the OMEGA grid is computed by First evaluating the 
top of the uppermost stretchable level (H, s ) above the lowest point in the domain (c.f Figure 21 right). 
This height is then used as a spherical top for the stretchable levels and the heights of each of the other 
vértices in the stretchable levels are computed by calculating a local geometric factor a, for each radial. 
Since H, s was computed using the lowest point in the terrain, these a¡ will always be less than CW- The 
resulting grid is terrain following at the surface, but has a constant mean sea level altitude at the top of 
the uppermost stretchable layer. 

An altemative setup method exists in which the user specifies the top of the total domain, and the 
OMEGA grid generator then calculates the geometric stretching factor. A final means of specification of 
the vertical grid structure is available to the more sophisticated user by editing the grid generator input 
file to provide the height of each vertical level above the surface. In this mode, the grid generator uses 
these heights to construct the radial through the lowest point on the surface grid, and then adjusts the 
other radiáis to produce a spherical top to the stretchable layers. 

Calculation of the vertical levels completes the construction of the OMEGA model grid. With the grid 
completed, all vital information is handed over to the model preprocessor. This ineludes the number of 
vértices, edges, and faces, and the arrays that contain the vertex locations, their connectivity, the edge 
connectivity, and the face connectivity. 
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5.3 Surface Characteristics Databases. 


The OMEGA Grid Generator accesses some surface datasets during the adaptation process (e.g. , terrain 
elevation and land/water fraction) while others are only interpolated after the grid has been defined (e.g., 
vegetation, soil moisture, sea surface temperature). Terrain and shoreline datasets have been described 
in Section 5.2.1. After the grid has been generated, the grid-generator reads the surface characteristics 
datasets and assigns valúes of sea surface temperature, landuse, sub-soil temperature, soil moisture, 
vegetation type and soil type to each grid cell. 

The use of an easily configurable adaptive grid requires a set of high-resolution surface databases, built 
in such a way that they can be efficiently accessed when a new grid is desired. OMEGA has eight major 
surface characteristic databases: soil type, climatological vegetation Índex, land use/land cover, 
climatological sea surface temperature, climatological subsurface temperature and climatological soil 
moisture. The first two characteristics are determined from global datasets, which have been widely 
used for numerical modeling. The others rely on newer datasets from various sources. For both 
vegetation index and land cover, two different datasets are provided. One dataset covers North America 
only at high resolution, and the other covers the entire globe, but with less resolution. Tableé lists the 
surface characteristics databases which are available for OMEGA. 

A one-degree (110 km) global soil type database was created from the Global Ecosystems Database 
(GED) CD-ROM (Webb et al., 1992). The source dataset consists of sepárate fields containing the 
amounts of clay, sand and silt in 12 soil types. The twelve types used in the OMEGA preprocessor are 
sand, loamy sand, silt loam, loam, sandy clay loam, silty clay loam, silt, clay loam, sandy clay, silty clay, 
and clay. 

Vegetation modifies both the thermal and moisture budgets of the Earth’s surface and thereby influences 
atmospheric circulations. In OMEGA, vegetation type is defined for each grid cell based on the Global 


Table 3. Datasets used to determine surface characteristics. 


Surface Property 

U.S. Data 

Global Data 

Soil Type 

Webb Global Soil Texture Class 
from Global Ecosystems Database (GED) CD-ROM 

1 degree resolution 

Climatological Vegetation 
Index 

USGS NDVI Annual Datasets 

2 minute resolution 

Land Use/ 

Land Cover 

USGS BATS Land Cover 

2 minute resolution 

Climatological Sea Surface 
Temperature 

Bi-weekly USGS SST 
Climatology 

12 minute resolution 

Monthly Global Vegetation 
Index (GVI) from GED 

10 minute resolution 

Climatological Subsurface 
Temperature 

Average Monthly Air 
Temperature from GED 

30 minute resolution 

Olson World Ecosystems 
BATS Land Cover from GED 
30 minute resolution 

Climatological Soil Moisture 

Budyko Soil Moisture Budget using 

Average Monthly Temperature and Precipitation from GED 

30 minute resolution 
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Vegetation índex (GVI). The GVI has a resolution of about 13 km. More recently, the USGS has 
published a 1 km Normalized Difference Vegetation Index (NDVI) over the contiguous U.S. Two 
climatological NDVI dataseis were created for OMEGA: a two minute (4 km) dataset covering the 
conterminous U.S., and a 30 minute (55 km) dataset covering the globe. A series of high resolution 
AVHRR dataseis (Eidenshink, 1992) over the conterminous U.S. has been published by the Earth 
Resources Observation System (EROS) Data Center of the U.S. Geological Survey. Dataseis beginning 
in 1990 are available, which contain gridded NDVI valúes, which have been composited over two 
weekperiods covering most of each year. The compositing process is an attempt to eliminate the effects 
of contamination by clouds. The 1990 through 1993 datasets were averaged to produce bi-weekly 
climatological datasets covering the conterminous U.S. with a resolution of 2 minutes (about 4 km). 
Global monthly NDVI datasets covering the period 1985-1990 were published on the Global Ecosystems 
CD-ROM, based on biweekly calibrated GVI data (EDC-NESDIS, 1992). 

These monthly datasets were averaged to form a set of monthly climatology files with a resolution of 
10 minutes (18 km). Similarly, two sepárate land use datasets were created. Each set uses the 
Biosphere-Atmosphere Transfer Scheme (BATS) set of 18 land cover categories, with an additional 19th 
category added for urban land cover. The categories are: crop/mixed farming, short grass, evergreen 
needleleaf tree, deciduous needleleaf tree, deciduous broadleaf tree, evergreen broadleaf tree, tall grass, 
desert, tundra, irrigated crop, semidesert, ice cap/glacier, bog or marsh, inland water, ocean, evergreen 
shrub, deciduous shrub, mixed woodland and urban. A 2 minute (4 km) land use dataset covering the 
U.S. was formed from a 1 km BATS dataset and AVHRR data. A 30 minute (55 km) global dataset was 
formed from the Olson (1992) World Ecosystems dataset on the Global Ecosystems Database CD-ROM. 
The 73 categories in the original Olson dataset were remapped to the 18+1 BATS categories. 

OMEGA uses a global sea surface temperature climatological dataset produced by the USGS 
(Schweitzer, 1993). The data has been processed into a bi-weekly global, 0.2 degree (about 20 km) 
dataset. This data does not extend beyond about 70 degrees latitude in either hemisphere. Some areas 
are not covered resulting in large lakes without data; for example the Caspian Sea in central Asia is 
completely missed. The original dataset was weekly; every other week was processed for OMEGA. 

A monthly global subsurface temperature climatology was built from monthly average surface air 
temperatures from the Global Ecosystems Database CD-ROM (Legates and Willmott, 1992). This data 
was built by using monthly average surface air temperature as an estimate of subsurface temperature. 

The dataset resolution is 30 minutes. Soil moisture is an important variable in numerical weather 
prediction but is difficult to determine. Budyko’s simple soil moisture budget as described by Sellers 
(1965) was used to generate soil moisture estimates from the temperature and precipitation data (Legates 
and Willmott, 1992). The monthly budget equation is given by: 

w 2 = Wj + r- E — S , (184) 

where W2 is the soil moisture valué at the end of a month, wi is the soil moisture at the beginning of the 
month, r is the precipitation, E is the evapotranspiration, and S is the moisture surplus term. The key 
term, evaporation, is based on a potential evaporation, a function of the monthly mean temperature. This 
simple approach is far from perfect, but seems to give reasonable results for a range of climatic types. 

An iterative application of the Budyko equations was used to obtain volumetric soil water contents that 
balanced the annual water cycle. Each of the surface characteristics databases are organized as sets of 
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latitude-longitude "tiles", which results in complete geographic coverage, efficient searching and logical 
filenames. The global-scale databases use 30 x 30 degree tiles, while the U.S. databases consist of 5 x 5 
degree tiles. FORTRAN logic was devised to allow databases at both resolutions (U.S. and global 
coverages) to be used simultaneously. The higher resolution data is applied to the OMEGA grid first, 
then the coarser resolution data is used in any grid cells which have not been filled with valúes. After all 
of the available databases have been searched, the database valúes which have accumulated in each 
OMEGA grid cell are averaged to obtain the final valué. In the case of the discrete fields (soil type and 
land use), the dominant (most frequent) valué in each cell is considered to represent the entire cell. If no 
database valúes are found for a grid cell, the cell is flagged and assigned a -99.0 valué. The PBL 
routines later assign valúes for these cells before the model run starts. Figure 22 shows some of the 
different data sets that are accessed by the OMEGA Grid Generator. 

5.4 Static Grid Adaptation. 

To the extent allowed by the minimum and máximum resolutions requested, the grid generator adapts 
the grid to features of the geography by creating small triangles in areas where gradients of the 
underlying data are high, while leaving triangles large in areas of small gradients. In the case of 
land/water boundaries, a "high" gradient is a non-zero gradient. 

The ratio of minimum to máximum resolution (roughly, máximum to minimum edge length) specified 
must be greater than two to enable adaptation to occur. If the minimum and máximum resolutions are 
closer than this, the grid generator will be unable to trisect cells without the resulting cells being too 
small. A minimum to máximum resolution of two results in a grid in which cell size throughout the grid 
appears uniform to the eye. If the minimum and máximum resolutions specified by the user differ by 
less than twenty percent, the grid generator will set the máximum resolution to eighty percent of the 
minimum resolution. This prevenís computationally costly "thrashing": repeated trisecting and 
recombining of cells in an attempt to satisfy resolution criteria that are too restrictive. 

The user can specify through the XOMEGA graphical user interface a point location around which 
refinement is required (c.f. Figure 23). The resolution and the radius of influence is also required as an 
input. XOMEGA generates a grid.ref file which is read by the grid-generator as an input for point 
refinement. The OMEGA grid-generator is configured to read in as many as fifty point sources. An 
example of grid.ref is given below: 

1 ! number of points 

10.25 35.25 40.00 1.00 ! longitude latitude radius(km) resolution (km) 

The grid-generator defines a Gaussian function around the point location using the input radius of 
influence. This function is used to provide weighting to the cell badness function. The cells cióse to the 
point source get larger weighting and are therefore tagged for refinement. The major advantage of point 
refinement is in reducing costly I/O time, while providing resolution exactly where it is needed. For 
large domains, using the point refinement option can generate high-quality grids relatively quickly. This 
technique can also be easily employed for refinement around line sources. Figure 23 shows the 
computational domain for the Middle East. The final grid has 4213 cells with edge lengths varying from 
5.33 km to 115.57 km. For this grid a single domain with nine point sources was defined. The grid.inp 
and the grid. ref files used to generate the grid are as follows: 
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(c) Land use. 


(d) Soil type. 



(e) Climatoligical vegetation 


(f) Climatological soil temperatura. 



(g) Climatological soil moisture. 


(h) Climatological sea surface temperatura. 


Figure 22. Datasets used to initialize OMEGA. 
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grid.inp: 

new 

oldgrid.dat 
grid.dgn 
35.00 
55.00 
22.00 
40.00 
40.00 
100.00 
5min 
5min 


new grid 

•k -k k -k k k k 

grid diagnostics file ñame 
minimum longitude (degrees); 
máximum longitude (degrees) 
minimum latitude (degrees) 
máximum latitude (degrees) 
minimum resolution (km) 
máximum resolution (km) 
terrain dataset resolution 
land/water dataset resolution 


smoothing = 


0 


! number 

of sub-domains 


30 


! total number of vertical levels 

28 


! number 

of stretched levels 

1 


! type of 

vertical grid 


30.00 


! height 

of the first vertical level (m) 

1.15 


! stretch 

ratio 


15 


! máximum 

gridgen iterations 

0 


! 0 = no 

; 1 = yes 

(create . srf. fi 

grid.ref: 

Q 





46.120 

28.330 

100.000 

7.000 

! Hafr-al-Batin 

50.150 

26.270 

100.000 

7.000 

! Dhahran 

51.570 

25.250 

100.000 

7.000 

! Doha 

43.520 

32.380 

100.000 

7.000 

! Ukhaydir 

46.720 

24.720 

100.000 

7.000 

! Riyadh 

39.700 

24.550 

100.000 

7.000 

! Medinah 

41.680 

27.430 

100.000 

7.000 

! Hail 

51.350 

35.680 

100.000 

7.000 

! Tehran 

46.380 

30.750 

100.000 

7.000 

! Khamisiyah 



Figure 23. OMEGA grid for the Middle-East región used for gulf war 
analysis. The grid has been refined around data sources. 
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Figure 24. Refinement to regional interim emissions inventory. 

Another example of the flexibility of the OMEGA grid is seen in Figure 24. This figure shows an 
OMEGA grid generated to provide higher resolution in the vicinity of emission sources. In this case, the 
grid adapted to a subset of the 1990 Regional Interim Emissions Inventory (sources producing more than 
25,000 tons per year (TPY) of S0 2 ). OMEGA, with dynamic adaptation, could be configured to provide 
high resolution not only in the source región, but also in those regions with high pollutant concentration. 
This is important in atmospheric chemistry situations where the lack of grid resolution can drastically 
change the solution due to the strong non-linearity of the physics This capability is discussed in detail in 
Section 9. 

5.5 Grid Editor. 

OMEGA Grid Editor has been developed to compliment the grid-generation process. Large amount of 
CPU time is required to read the terrain and shoreline datasets. To avoid going through the whole 
process in order to make minor changes to the existing grid, which is an inefficient use of resources, 
SAIC developed the Grid Editor to perform a variety of functions on an existing grid. Using the Grid 
Editor, areas can be defined for refinement or coarsening, the entire grid can be regenerated with 
different resolution without creating the terrain and shoreline datafiles from scratch, and individual cells 
can be manipulated by deleting vértices or bisecting boundary edges. Bisection of boundary edges is 
especially useful, since the chances of creating bad cells on the boundary are high. In addition, overall 
grid-quality can be improved by relaxing the vértices. The Grid Editor is a menú driven program which 
should be used interactively with XGRID. The main menú of the Grid Editor is shown below: 
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* OMEGA GRID EDITOR, VERSION 3.5 

* * MAIN * * 


r -> 

Refíne Grid. 

c -> 

Coarsen Grid. 

b -> 

Bisect Boundary Edge. 

v -> 

Delete Vertex. 

X -> 

Relax Grid. 

s -> 

Change Resolution. 

u -> 

Undo. 

q -> 

Quit. 


5.6 Dynamic Grid Adaptation. 

The grid adaptation described so far is “static” adaptation, implemented in the grid generator prior to the 
beginning of an OMEGA run. The grid in the “static” adaptation case is adapted a prioi to resolve static 
features such as terrain gradients and land-water boundaries, or any other feature that the user ineludes in 
the adaptation scheme. The grid does not change during the course of the OMEGA run. The latest 
versión of OMEGA also has dynamic adaptation capabilities. This new feature is discussed briefly here 
(for details see Section 9). 

OMEGA has the ability to adapt its grid during a simulation to different criteria, such as, frontal activity 
and convection or a pollutant plume. The underlying philosophy of dynamic adaptation is to provide a 
refined grid to resolve important physical events and features as the simulation is in progress. The 
dynamic adaptation consists of four major steps: 

1) At a predetermined time step specific variables or their gradients are evaluated to see if they meet 
the adaptivity criteria. 

2) The mesh is refined where these criteria are satisfied. 

3) The physical variables are interpolated to new cell centers, and 

4) The mesh is coarsened where the criteria are not met. 

This process can be accomplished in several ways. For example, if there are sharp discontinuities such 
as propagating fronts the mesh is refined ahead of the front so that a refined mesh already exists when 
the shock arrives at a specific location. As the front propagates, the refined mesh at the previous 
location of the front is coarsened to the resolution specific to the background mesh. This type of 
coarsening and refining is performed at predetermined time-step increments and is useful for tracking 
frontal or shock activity. The simulation is allowed to advance a number of timesteps defined through 
an input file without adaptation and then refine-coarsen eyele is performed. However, if we choose to 
dynamically adapt the mesh to resolve spontaneous events we need to check predetermined physical 
variables at prescribed time intervals to see if they exceed a threshold or satisfy a criterion. We then 
perform a refine-coarsen eyele. After the grid has been modified all the physical variables are 
interpolated to the new mesh and the diagnostic variables are calculated. 
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(a) Original grid. 


(b) Higher resolution over NW región. 



(c) Higher resolution over NE región. (d) Higher resolution over SE región. 


Figure 25. Test of the dynamic adaptation routine. The adaptation was to the 
presence of test partióles (shown in red) initialized in the Southwest 
comer of the grid and traveling on a prescribed trajectory. 

Figure 25 shows a simple test of the dynamic adaptation routines. The domain was centered on the Four 
Comers area. The adaptation criteria were set to point locations. A partióle traveling on a prescribed 
trajectory was initialized in the Southwest comer of the domain. Only the refinement cycles were 
invoked during this mn. 





Section 6 

Initial and Boundary Conditions 


The function of the data preprocessor is to convert a diverse mixture of real-time atmospheric and 
surface data for input into an OMEGA simulation. This process produces two dataseis. One dataset 
specifies the initial valúes of all the model’s dependent variables as well as the surface-based model 
parameters. A second dataset specifies the time dependent behavior of the model’s dependent variables 
along the lateral boundaries of the simulation domain. The general structure of the OMEGA 
preprocessor is depicted in Figure 26. 

The first module is INGEST. Its function is to ingest atmospheric data (e.g. surface or rawinsonde 
observations) from a real-time data stream and distribute it to a series of files (RTData) organized by 
type of data and time. INGEST continuously executes on the Computer system that is designated for 
data ingestión and is not part of the model run stream. 

Execution of the OMEGA grid generator is the first step in a preprocessor run. The latest versión of 
GRIDGEN also ingests data from the surface characteristics database and generates a specification of all 
surface properties for each grid cell. This capability is now operational in the latest versión of OMEGA. 

The DASH module reads real time surface data from both hourly reporting surface stations and daily 
reporting hydrological stations to maintain an updated global database of snow cover, soil temperature, 
and soil moisture data. DASH is not currently being used as part of the operational OMEGA system. 

The PREPDAT module selects data stored by the INGEST and DASH modules based upon the user- 
specified date and time. PREPDAT then error-checks the data and writes it out to a set of files written 
in standard OMEGA input format. This ineludes a "GUESS" file that contains a 3-dimensional first- 
guess of all dependent atmospheric variables interpolated from forecast model output. Other files 
inelude: 1) observational data files ("SFC", "RAOB" and "SST"); and 2) hydrological data files 
("TSUB", "WSHL", "WDEP" and "SNGR"). 

GENINIT assembles the atmospheric, and time varying surface property data into a single initialization 
file called "INIT". GENINIT also writes a "BC" file for the 3-D atmospheric variables at the initial 
time. 

ANALYSIS accepts processed observational data from PREPDAT files and performs an objective 
analysis of the data to update valúes of the dependent atmospheric variables at the centroids of all grid 
cells. ANALYSIS then outputs a modified "INIT" file. 

Finally, the preprocessor ingests forecast data from another atmospheric model and interpolates it to the 
lateral boundary cells of the OMEGA domain. This function is controlled by RUNPREP and is 
accomplished by running PREPDAT and GENINIT in Boundary condition mode once for each 
boundary condition time. GENINIT vertically interpolates the guess generated by PREPDAT to the 
OMEGA grid. The result is one "BC" file valid at a given forecast time for each time that PREPDAT 
and GENINIT are executed. 
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OMEGA receives from the preprocessor initial and boundary valúes for the latitudinal, longitudinal and 
vertical components of the wind velocity, pressure, potential temperature, density and mixing ratios of 
water vapor, cloud water, rain water, cloud ice and snow. In addition, OMEGA receives the initial 
valúes of skin and subsoil temperature; shallow and deep soil moisture; snow cover; and vegetation 
Índex. The input module of OMEGA converts the latitude-longitude-altitude coordinates into Cartesian 
coordinates. 

The files "BC.tag#” and "INIT.tag” are produced by the preprocessor, where tag represents the unique 
preprocessor run identifier. "BC.tag#" contains boundary condition information at time"#" starting with 
"0" at the initial time. "INIT.tag" contains the initial conditions at the OMEGA grid cells. 

6.1 Data Ingestión. 

Real-time gridded and observed atmospheric data is currently acquired from the Fleet Numerical 
Meteorological and Oceanographic Center (FNMOC). Until recently, gridded and observed data sets 
were also acquired from the National Centers for Environmental Prediction (NCEP) using their Internet 
ftp site. 

6.1.1 NCEP Internet Data Feed. 

NCEP stores their numerical model data in the World Meteorological Organization (WMO) Code FM 
92 GRIB format (Ed. 1), and their observational data in the WMO Code Form FM 94 BUFR format 
(Ed. 2). Observational data that are currently stored in NCEP Office Note formats are gradually being 
converted to BUFR format. Ultimately, all of the observational data will be stored in BUFR format. 

The packed binary BUFR format is more efficient in disk space utilization than the intemal NCEP 
formats. This will impro ve transfer times as well. 

The GRIB messages are comprised of packed binary data which is divided into six different sections as 
follows, some of which are optional: 

Sec. 0 Indicator Section 

Sec. 1 Product Definition Section (PDS) 

Sec. 2 Grid Description Section (GDS) - optional 
Sec. 3 Bit Map Section (BMS) - optional 
Sec. 4 Binary Data Section (BDS) 

Sec. 5 HIT (ASCII) 

The Indicator Section begins with the ASCII string GRIB’ which defines the beginning of the GRIB 
message, followed by the total length of the GRIB message and the Edition number of the GRIB 
message. The next section is the Product Definition Section that is usually 28 octets long, but may be 
longer. It contains information such as the numerical model that created the data, the center 
identification, grid identification, and date and time of forecast. It also contains flags specifying whether 
the next two sections are present. Section 2, the Grid Description Section, is optional and contains a 
grid description for grids not already defined as NCEP storage grids. Section 3, the Bit Map Section is 
also optional, and contains either a bit map for the data or a reference to a pre-defined bit map. Section 
4 the Binary Data Section, contains the binary scaling information which is used to reconstruct the 
original model data from the packed binary data in the GRIB file. The last section is four octets long 
and defines the end of the GRIB message with the ASCII string HIT. For further information on the 
GRIB message format, refer to GRIB Ed. 1 (Office Note 388), (Stackpole, 1994). 
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The BUFR file format is also a packed binary format, but is less structured than the GRIB format, and is 
more self-defining. Since there is a wide variety of data which may be stored in the BUFR format, from 
surface observations to aircraft and ship observations to radar and profiler data, each field of data stored 
in the BUFR format is stored with an accompanying numerical field called a descriptor. The descriptor 
is matched up against a table (TABLE B) which provides a set of information about that field, including 
the variable ñame, units, reference and scaling information and the number of storage bits required for 
that field. The BUFR decoder uses this information to correctly decode the data. 

Decoders were developed for both the GRIB and BUFR data files, using a combination of c-code to read 
in the binary data, and FORTRAN modules provided by NCEP to decode the binary data records and 
retum arrays filled with the unpacked data, which are then written out to ASCII files. The GRIB 
decoder uses an option file configured by the user to process the raw data file and save only the model 
grids and variables requested by the user. Often GRIB files contain data from more than one model and 
for many different variables. This flexibility to save only the data needed results in smaller files since 
decoding and saving entire GRIB messages can result in very large files. The BUFR file decoder 
doesn’t use an option file, since these files are usually smaller, and often contain just one data type at one 
observation time. Table 6.1 contains a list of NCEP files which the GRIB and BUFR decoders can 
process for input into the preprocessor. For the GRIB files, MM refers to the initialization time of the 
model simulation and is usually either 00 or 12. NN and NNN are the variable forecast times of the 
NCEP model runs relative to the model initialization time. LL refers to the three-hourly observation 
time of the data and ranges from 00 to 21 in three-hour increments. YYMMDDHH refers to an eight- 
digit date made up of the year (YY), month(MM), day (DD) and hour (HH). A represents a single 
character abbreviation for a variable: "R" for relative humidity, "U" for u wind component, "V" for v 
wind component, "H" for geopotential heights or "T" for temperature. 

6.1.2 FNMOC Data Feed. 

There are several global and regional model datasets available in real-time from FNMOC. There are 
two NOGAPS global model grids, which OMEGA is capable of using for initial conditions, and there 
are also several regional COAMPS model grids available, which will be used in future versions of 
OMEGA. They also provide rawinsonde data in the standard TTAA/TTBB format and surface 
observations in the METAR format. Table 4 contains a list of files available from the FNMOC data 
feed. 

For the GRIB files, X refers to a single letter that defines the model pressure level, W refers to a two- 
digit number that defines the variable type, and MMM refers to the variable forecast times of the 
NOGAPS model runs relative to the model initialization time. YYMMDDHH refers to an eight-digit 
date made up of the year (YY), month(MM), day (DD) and hour (HH). YYYYMMDDHH refers to a 
similar ten-digit date, but the year ineludes the century (eg. 1997). A represents a single character 
abbreviation for a variable: "R" for relative humidity, "U" for u wind component, "V" for v wind 
component, "H" for geopotential heights or "T" for temperature. 

6.1.3 Air Forcé Weather Agency (AFWA) Data Feed. 

An altemative access to real time atmospheric data for DSWA OMEGA applications is through the 
databases at the AFWA in Omaha, Nebraska. It is anticipated that most of the data Usted in Table 4 will 
also be able to be obtained through AFWA. The capability to decode and utilize the AFWA data is not 
included in versión 3.2.6 of the OMEGA preprocessor. It will be incorporated into subsequent versions. 
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Table 4. Summary of NCEP GRIB and BUFR data files ingested by PREPDAT. 


GRIB 

or 

BUFR 

Data Type 

Input File Ñame 
(stored at NCEP) 

Output File Ñame 
(read by PREPDAT) 


Gridded Data 

GRIB 

MRF 

CED1 .FNN.TMM.MRF 

YYMMDDHH*.mrfA 

GRIB 

NGM 190 km 

CED1 .RFNN.TMMZ.RGL 

YYMMDDHH*.ngmB 

GRIB 

NGM 90 km 

CED 1 .RFNNXP.TMMZ.RGL 

YYMMDDHH*.ngmG 

GRIB 

ETA 190 km 

CED 1 .LFNN.TMMZ.ETA 

YYMMDDHH*.etaB 

GRIB 

ETA 90 km 

CED 1 .LFNN.TMMZ.ETA 

YYMMDDHH*.etaG 


Rawinsonde Data 

BUFR 

Rawinsonde 

SNDDATA.BUFR.T00Z.MRF 

YYMMDDHH.snd 


Surface Data 

BUFR 


SFCSHP.NCBUFR.LATEST(CLLZ) 

YYMMDDHH.sfcs 

BUFR 


SFCSHP.NCBUFR.LATEST(CLLZ) 

YYMMDDHH.sstb 


Table 5. Summary of FNMOC data files ingested by PREPDAT. 


Data 

Format 

Data Type 

Input File Ñame 
(from raw FNMOC feed) 

Output File Ñame 
(read by PREPDAT) 


Gridded Data 

GRIB 

NOGAPS 2.5 deg 

XVY $MMMGRB 

YYMMDDHH*.fncH 

GRIB 

NOGAPS 1.0 deg 

XW1MMMGRB 

YYMMDDHH*.fncX 


Rawinsonde Data 

TTAA 

Rawinsonde 

YYYYMMDDHHZ.wmo 

YYMMDDHH.raob 


Surface Data 

METAR 

Land observations 

YYYYMMDDHHZ.sao 

YYMMDDHH.sao 


6.1.4 Other Data Types. 

Gridded and observational data types can be read in several formats that are archived by various 
govemment agencies including NCAR and NCEP (detailed in Section 6.4 describing the PREPDAT 
module); however some data types deserve special mention. The OMEGA preprocessor can read 
gridded output from previous OMEGA simulations or it can create a first-guess from a single tabular 
format hand-entered sounding when gridded data is unavailable. 

6.2 Processing of Surface Characteristics Databases. 

With the exception of the land/water and terrain elevation databases each of the other surface 
characteristics databases are organized as sets of latitude-longitude "tiles", which results in complete 
geographic coverage, efficient searching and logical filenames. The global-scale databases use 30 x 30 
degree tiles, while the U.S. databases consist of 5 x 5 degree tiles. 
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FORTRAN logic was devised to allow databases at both resolutions (U.S. and global coverages) to be 
used simultaneously. The higher resolution data is applied to the OMEGA grid first, then the coarser 
resolution data is used in any grid cells which have not been filled with valúes. After all of the available 
databases have been searched, the database valúes which have accumulated in each OMEGA grid cell 
are averaged to obtain the final valué. In the case of the discrete fields (land/water, soil type and land 
use), the dominant (most frequent) valué in each cell is considered to represent the entire cell. If no 
database valúes are found for a grid cell, the cell is flagged and assigned a "last resort" default valué 
later in the preprocessor sequence. 

6.3 Processing of Real-Time Hydrology Data. 

The DASH (Daily Analysis of Soil and Hydrology) program can be used to process real-time reports of 
precipitation, temperature, and snow cover in order to produce analyzed fields of model-ready 
parameters. Soil moisture at various depths, subsurface temperature and snow cover are output in tiled 
latitude-longitude files. 

DASH reads hydrological data in four specific formats. DASH performs objective analyses of daily 
hydrological variables from a set of sites and produces a set of gridded latitude-longitude tiles. These 
primary variables may be processed to obtain soil moisture, snow cover and subsurface temperature, 
which are needed for input into OMEGA. Error checking is also performed by DASH. A simple 
inverse-distance weighting scheme is used for the objective analyses. 

Daily precipitation is analyzed and then used by the Antecedent Precipitation Index (API) scheme to 
estímate soil moisture. The retention coefficient used in the API scheme varíes as a function of soil 
depth. The subsurface temperature is estimated by averaging the surface air temperature over a three 
day period. Snow cover is calculated by analyzing each day’s snow cover observations. An API-type 
algorithm is also applied to gradually deplete the previous day’s snow cover. 

The primary variables are: 

(1) máximum surface air temperature, 

(2) mínimum surface air temperature, 

(3) total daily precipitation, 

(4) máximum soil temperature, 

(5) mínimum soil temperature, 

(6) depth of soil temperatures, 

(7) soil moisture, 

(8) depth of soil moisture, 

(9) evaporation, 

(10) solar radiation; 

(11) duration of vegetation wetness, 

(12) total snow on the ground, 

(13) new snowfall, and 

(14) average daily temperature. 

The secondary variables are: 

(1) Antecedent Precipitation Index (API) at 5 cm and 30 cm, 

(2) estimated volumetric soil moisture contents, 

(3) average surface air temperature, 

(4) 24-hour precipitation, and 

(5) snow cover. 
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6.4 PREPDAT Module. 


The PREPDAT module of the preprocessor serves as the interface between the many forms of 
atmospheric data which can be utilized by the OMEGA initialization process and the modules which 
actually use the data to create a three-dimensional initialization dataset for OMEGA. PREPDAT is 
structured to ingest data from the RTData set of files generated by the INGEST module as well as other 
archived atmospheric data formats. This permits the OMEGA initialization process to be executed for 
both research and real-time simulations. After ingesting the data PREPDAT performs a series of simple 
error-checking procedures. More stringent quality control is performed by ANALYSIS. Once the error- 
checking is complete each type of data is written to a file in a standard format. 

PREPDAT ingests rawinsonde, surface and sea surface temperature observations, first guess gridded 
data and gridded hydrological data. These data sources are described in Table 6. 

6.4.1 Gridded First-guess Data. 

Currently, OMEGA can ingest the following types of first-guess data: 

1. Data from several archive databases maintained by various govemment agencies including 
LFM II/Bureau of Reclamation format, NGM data in UNIDATA NetCDF format. Global 
Optimum Interpolation from NCAR, data on the ETEX experiment grid from ECMWF, and 
ECMWF 2.5 degree lat/lon global analysis from NCAR; 

2. Data decoded from NMC’s GRIdded Binary (GRIB) format available in real time by FTP 
from NMC including MRF output, and NGM and ETA output on several different grids; 

3. Data from FNMOC’s models also decoded from GRIB format; 

4. Data from AFWA including HIRAS and HYSPLIT; 

5. Data from a previous OMEGA simulation; and 

6. A single hand entered sounding. 

The user must choose exactly one gridded data source from the list in the user options file 
"prepomega.opt". PREPDAT processes the gridded data in the following manner: 

1. The data file is read; 

2. The data are interpolated horizontally to grid cell centroids and vertically to the pressure 
levels set by the user in "prepomega.opt"; 

3. The missing valúes and levels are filled in; and 

4. The data are written to the "GUESS" file. 

6.4.2 Rawinsonde Data. 

Rawinsonde data is optional but recommended for every OMEGA initialization. The types of 
rawinsonde data that can be ingested by OMEGA are listed in Table 6. The user can choose one or 
more rawinsonde data source from the list in the user options file "prepomega.opt". PREPDAT 
processes the rawinsonde data in the following manner: 

1. The data file is read; 

2. At each rawinsonde site, the mandatory level, significant level and wind level data are 
merged into a single sounding; 

3. Each sounding is checked for errors such as temperature spikes, superadiabatic layers and 
spikes in the wind hodograph, the errors are corrected, if possible, or the data is eliminated; 
and 

4. The data are written to the "RAOB" file. 
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Table 6. Summary of data types which can be ingested by PREPDAT. 



Data Type 


Gridded Data 


OMEGA output 


LFM II/Bureau of Reclamation format 


NGM (UNIDATA NetCDF format)_ 


Global Optimum Interpolation 


ECMWF 2.5 degree lat/lon global analysis (for validation cases). 


Decoded GRIB data from NCEP - MRF 


Decoded GRIB data from NCEP - NGM or ETA 190 km 


Decoded GRIB data from NCEP - NGM or ETA 90 km 


Decoded GRIB data from NCEP - ETA 81 km Lambert Conf. (Grid L 
or 211) _ 


Decoded GRIB data from FNMOC - NOGAPS 2.5° global grid 


Decoded GRIB data from FNMOC - NOGAPS 1.0° global grid 


ETEX Experiment grid from ECMWF 


Uniform guess from hand-entered soundin 


HIRAS data from AFWA 


HYSPLIT data from AFWA 


Rawinsonde Data 


"<caseid>DHHMMz.out" 


"GborYYMMDDHH" 


"MM * DDHH.cdf" 


"GgoiYYMMDDHH" 


"GecmYYMMDDHH" 


"YYMMDDHH*.mrfA" 


"YYMMDDHH*.ngmB" 

"YYMMDDHH*.etaB" 


"YYMMDDHH*.ngmG" 

YYMMDDHH*.etaG" 


"YYMMDDHH*.etaL" 


"YYMMDDHH.fnlH" 


"YYMMDDHH.fnlX" 


not implemented 


"RAOBFG" 


"GhirYYMMDDHH" 


user specified in option file 


Decoded from NCEP Office note 29 real-time format 

"RnmcYYMMDDHH" 

Decoded from NCEP Office note 29 archive format (for validation 
cases). 

"RnmcYYMMDDHH" 

Decoded from NCEP BUFR format 

"YYMMDDHH.snd" 

TTAA/TTBB/PPBB format from NCEP or FNMOC 

"YYMMDDHH.raob" 

DSWA historical format 

"YYMMDDHH.rhst" 

Air forcé archive format 

"RafaYYMMDDHH" 


Surface Data 


Decoded From NCEP Archive format (for validation cases). 


SAO data 


DSWA historical format 


Land observations decoded from NCEP BUFR format 


Office note 29 land surface observations (real-time format). 


Office note 29 land surface observations (archive format). 


METAR format 


Buoys and Ship Reports 


Ship observations decoded from NCEP BUFR format 


Sea Surface Temperature Data 


Office note 29 land surface observations (real-time format). 


Weekly OI analysis from NCEP 


"SmasYYMMDDHH" and 
"YYMMDDHH.SEA" 


"YYMMDDHH.sao” 


"YYMMDDHH.shst" 


"YYMMDDHH.sfcl" 


" Y YMMDDHH .sfcl" 


YYMMDDHH.sfcl" 


"YYMMDDHH.sao" 


"YYMMDDHH.SEA" 


"YYMMDDHH.sfcs" 


"YYMMDDHH.sfcl" 


"oi .mean .bias. YYMMDD" 

























































6.4.3 Surface Data. 

Surface data is optional but recommended for every OMEGA initialization. The types of surface data 
which can be ingested by OMEGA are listed in Table 6. The user can choose one or more surface data 
source from the list in the user options file "prepomega.opt". PREPDAT processes the surface data in 
the following manner: 

1. The data file is read; 

2. Each observation is checked for obvious errors such as out of range valúes with erroneous 
data eliminated; and 

3. The data are written in a standard tabular format to the "SFC" file. 

6.4.4 Sea Surface Temperature Data. 

Sea surface temperature (SST) data is optional but recommended for every OMEGA initialization if the 
domain ineludes large water bodies. The types of SST data which can be ingested by OMEGA are listed 
in Table 6. The SST can be initialized using a biweekly climatological database or a new weekly real- 
time SST analysis which recently became available on the NCEP NIC server. This optimum 
interpolation (OI) analysis is produced on a one-degree grid using in situ and satellite SST’s plus SST’s 
simulated by sea-ice cover. This dataset is usually available by anonymous ftp each Monday for the 
preceding 7-day period. 

The user can choose one or more SST data source from the list in the user options file "prepomega.opt". 
PREPDAT processes the SST data in the following manner: 

1. The data file is read; 

2. Each observation is checked for obvious errors such as out of range valúes with erroneous 
data eliminated; and 

3. The data are written in a standard tabular format to the "SST" file. 

6.4.5 Hydrological Data. 

Hydrological data is optional but recommended for every OMEGA initialization. Climatological valúes 
of subsurface temperature and soil moisture are available in the on-line databases. GENINIT 
automatically uses the climatological data wherever there is no real-time data available. If no real-time 
data is available, climatological data is used everywhere. The types of hydrological data which can be 
ingested by OMEGA are listed in Table 6. The user can choose one or more hydrological data sources 
from the list in the user options file "prepomega.opt". PREPDAT processes the hydrological data in the 
following manner: 

1. The appropriate real-time data files are read; 

2. The data are averaged and/or interpolated to the OMEGA grid - In regions where the 
OMEGA grid is coarse relative to the data, averaging is performed / if the data is coarse 
relative to the OMEGA grid, interpolation is performed and 

3. The data are written out to the files: TSUB (subsurface temperature), WSHL (shallow soil 
moisture), WDEP (deep soil moisture) and SNGR (snow cover). 

6.5 GENINIT Module. 

The GENINIT module of the data preprocessor assembles the atmospheric and surface property data. 
GENINIT reads the SSTclim, NDVIclim, WSOILclim and TSUBclim files. It then reads the GUESS, 
TSUB, WSHL, WDEP and SNGR files from PREPDAT. GENINIT writes out an "INIT" file that 
specifies both atmospheric and surface property variables and a "BC" file that specifies only the 
atmospheric variables at boundary cells. 
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6.5.1 Surface Property Variables. 

GENINIT assigns the valué of the hydrology variables at each grid cell in the following order of 
preference: 

1. Non-flag real-time valué from the DASH generated files: TSUB, WSHL, WDEP or SNGR; 

2. Climatological valúes of soil moisture and subsurface temperature are used when the DASH 
valúes are unavailable - valúes from the WSOILclim file are used for both shallow and deep 
soil moisture; and 

3. If neither DASH ñor climatological valúes are available, the default specified by the user in 
"prepomega.opt" is used. 

NDVI and SST are assigned a default valué if the climatological valúes are unavailable. SST can be 
updated with observations in ANALYSIS, if they are available. 

6.5.2 Atmospheric Variables. 

The atmospheric variables are vertically interpolated from the constant pressure surfaces provided by 
PREPDAT to the centroids of the cells specified in the "grid.out" file. The variables are interpolated 
linearly with respect to height. 

6.6 ANALYSIS Module. 

Before a numerical simulation can be executed the many forms of raw initialization data must be 
combined into a single dynamically consistent description of the State of the atmosphere at the time of 
initialization. This involves: (1) filtering erroneous or inconsistent data from the initialization 
procedure; (2) obtaining initial valúes of atmospheric and surface variables at model grid point locations 
from an array of often irregularly-spaced data from different observing Systems (objective analysis); and 
(3) adjusting the initial State to remove unrealistic relationships among the variables caused by 
erroneous, unrepresentative or sparse data (initialization). It is the function of the ANALYSIS module 
of the preprocessor to perform these functions and update a model initialization file "INIT". 

The core of the OMEGA analysis and initialization procedure consists of seven steps: 

1) the ingestión of the available data prepared by the PREPDAT module and the preliminary 
"INIT" file created by the GENINIT module; 

2) the objective analysis of sea surface temperature data; 

3) the processing and quality control of surface and upper air data; 

4) the objective analysis of surface and upper air data; 

5) adjustments to insure that the analyzed data meets desired constraints; 

6) the generation of objective analysis performance statistics; and 

7) the output of the initialization data to a modified "INIT" file. 

6.6.1 Setting Up the Grid and Reading the Data. 

The first step in the analysis and initialization procedure is the input of user options from the 
"prepomega.opt" file. These options permit the user to select many of the attributes of the analysis and 
initialization procedure: (1) the choice of data that is included in the analysis procedure; (2) the degree 
of geostrophic coupling between the mass and momentum fields in the multivariate optimal 
interpolation (OI) scheme; and (3) the e-folding radius of the correlation functions used in the OI 
scheme. The next step in the procedure is the input of the data from the files that define the grid that 
will be used for the model simulation. The information that is ingested ineludes the latitude, longitude 
and elevation of each cell vertex, the terrain elevation, the land and water designation of each grid point, 
the land use classification, the climatological NDVI and the climatological sea surface temperature 
(SST). 
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After the input of the grid description and user option data, a sequence of subroutines which ingest 
rawinsonde, surface and sea surface temperature data from PREPDAT are executed. The flags set by 
the specifications in the "prepomega.opt" file determine which of these subroutines are executed and, as 
a result, which dataseis are ingested into the ANALYSIS module. The preliminary "INIT" file from 
GENINIT is also read at this time. 

6.6.2 Sea Surface Temperature Analysis. 

Once all the data has been ingested, the next step is to analyze the SST. The SST first-guess is prepared 
in the following manner: If a first-guess is available, it is "smeared" onto land grid cells that are within 
2 grid cells of a water point. This improves the calculation of differences between the first-guess and 
the observations. If no SST first-guess is available, the first-guess level one air temperature is used. 

The SST guess data is a static database of biweekly 30 year average sea surface temperature data 
interpolated to the model grid. Current SST data is generally considered to be observations that were 
made up to 7 days prior to the time of initialization. The use of this time period is based on the premise 
that sea surface temperatures change relatively slowly and that SST reports from many ships are at 
infrequent intervals. The analysis consists of a two-dimensional versión of the OI scheme. 

6.6.3 Observational data processing and quality control. 

The next step is to process the rawinsonde and surface data. Each rawinsonde sounding is searched for 
gaps that are larger than 50 mb. These gaps are filled through interpolation. Next, the first-guess is 
interpolated to the location of each rawinsonde observation. This valué is then subtracted from the 
observation to produce a residual or "observation increment". Finally, the rawinsonde observation 
increments are interpolated to mandatory pressure levels and the root mean square observation 
increment is calculated for winds and heights at each level. These valúes are used by the multivariate OI 
scheme as well as in the verification of the analysis. This is done by objectively analyzing not only at 
the model grid points but also at the rawinsonde sites at mandatory levels and at the surface observation 
sites. The RMS differences between the analysis and the observations is then compared to the RMS 
differences between the first-guess and the observations. 

Surface data are converted to observation increments as are rawinsonde data. A search is then made for 
duplícate observations which are flagged. Then the surface observation list is Consolidated and flagged 
observations are removed. 

At this point quality control is done for both surface and rawinsonde data. The quality control algorithm 
is essentially equivalent to the horizontal check component of Gandin’s (1988) complex quality control 
algorithm. The 3-D OI scheme is used to analyze the observation increments O-B (where is O is the 
observed valué and B is the first-guess interpolated valué) from nearby observations not including the 
observation itself to the location of each observation. The observation increment at each site is then 
subtracted from the analyzed valué to determine how well each observation increment agrees with its 
neighbors. This difference between the observation increment and the valué of nearby observation 
increments objectively analyzed to the observation location is called D. If the valué of D for a given 
observation is larger than a threshold, the observation is thrown out. The current error-thresholds are 
shown in Table 7. Up to 8 observations are allowed in each matrix for the rawinsonde QC algorithm 
and 6 for the surface QC algorithm. QC for rawinsondes ineludes height, relative humidity and 
horizontal wind. In addition, the level 1 height correction is checked against height corrections for the 
complete set of surface and rawinsonde observations. If the level 1 height corrections differ from the 
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Table 7. 


Observation Increment Error Thresholds. 


Variable 

Surface 

Rawinsonde 

Height (m) 

20.0 

20.0 at surface, increasing linearly to 50.0 
at top of model domain. 

Relative humidity (fraction) 

0.3 

0.3 

Horizontal wind vector 
(% of wind speed) 

50.0 % or 5 ms -1 , 
whichever is greater 

22.0% 

Temperature (K) 

4.0 

4.0 


analyzed valué by more than the threshold, all height reports from the sounding are tossed. If the 
disagreement is greater than 200 m, the entire sounding is tossed. For surface data, horizontal wind 
vector, height, temperature and relative humidity corrections are subject to quality control. 


After the quality control is performed, the last step before the multivariate OI is the combination of 
surface observations into "superobservations" to help improve the conditioning of the correlation 
matrices which are solved by the 3D-OI scheme. The number of surface observations within 45 km of 
each observation are counted. Those with the most nearby observations are processed first. A 
simplified versión of the OI scheme is used to interpólate the observation increment for each surface 
variable to the site of the central observation. This valué is taken as the superobservation and the others 
are eliminated. This procedure is iterated until no two surface observations are within 45 km of each 
other. 


6.6.4 Three Dimensional Multivariate Optimal Interpolaron Scheme. 

The basic principie in the OI scheme is that an estímate of the valué of a variable Ai at a set of grid 

points (i = 1,2,.I) can be created from a set of "J" meteorological observations Oj (where j=l, 2, ...,J) 

distributed irregularly throughout the analysis domain by forming an "optimal" linear combination of the 
"K" observations closest to the grid point. Following Daley (1991), this approach yields the following 
equation for the analyzed valué of Ai at the grid point: 


A 


=*¡ + X w * 


k =1 


rms, r „ -i 
rms k 


(186) 


where Ai = analyzed valué at model grid points 

Bi = first guess or background valué at model grid points 

0^ = meteorological observations 

K = total number of observations allowed to influence grid point i, 

Wk = weight associated with observation k 

rmsi = root mean square first-guess error for variables of the same type as the analysis 
variable 

rmsk = root mean square observational error for variables of same type as observation k, 

Bk = first-guess or background valué at observation locations 
Ok = observed valué. 

The rms errors are computed by averaging horizontally over many events. In the OMEGA preprocessor, 
both of the rms errors are approximated by the rms difference between the observation and the first- 
guess at a given pressure level. In the univariate case, they cancel each other, since the analysis variable 
is always the same as the observation variables. The method to obtain the "optimal" linear weighting 
coefficients for each of the observations is based on Gandin’s (1963) theory of optimum interpolation. 
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The weighting coefficients are determined by attempting to minimize the root mean square error 
(T k - of the interpolated valúes (i.e. the analyzed valúes) where T represents the true valué not 

including signáis at scales that cannot be resolved by the observing network and the overbar denotes an 
ensemble average of many cases. The mínimum of the root mean square error of the analyzed valúes 
can be found by taking a pardal derivative of Equation 6.2 with respect to the weighting coefficients, 
and setting it to 0 to form the set of equations: 
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( 188 ) 


where K = total number of observations allowed to influence grid point i; 

T = true valué not including signáis at scales that cannot be resolved by the observing network; 
B-T = first-guess error; 

O-T = observational error; 

P {b-t) - correlation between the first-guess errors at observing locations K and L; 

P [o-t) ~ correlation between the observational errors at observing locations K and L; 

P (b-t) =correlation between the first-guess errors at observing location K and grid point i; 

e = ratio of the mean observational error to the mean first-guess error, assumed to be invariant 
in space. 


The observational error ineludes two parts, the instrument error and the error of representativeness 
produced when information at very small wavelengths is aliased to the larger wavelengths due to the 
coarseness of the observing network. Essentially the analyzed valué is the sum of the first-guess and the 
sum of the weights multiplied by the difference between the observation and the first-guess at each 
observation location. The weights are determined by solving the system of equations represented by the 
matrix problem in Equation 188. OI is superior to successive correction schemes such as the Bames 
(1964) scheme since it takes into account the locations of observations relative to each other and 
effectively reduces the weight given two nearly collocated observations. It also takes into account the 
relative reliability of the observations relative to the first-guess. 

The two dimensional matrix in Equation 188 can be simplified in several ways. First, notice that it is 
symmetric about the diagonal since P( B -t) u + e P(o-T) k , ~ P(BT) lk + e P(o-r) lk • Also, with a few exceptions, 

P(o-T)¡ k ~ 0 when j ^ ^ and P( 0 -T) ¡k = 1 w hen l=k, since most observations are independent of each other. 

The two most notable exceptions are observations at different levels in the same rawinsonde sounding 
and observations from the same satellite-based sensor. Since P( B - T ) kt = 1 when l=k, the diagonal terms 

in the matrix simplify to 1+e. 
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The true valué, T, is unknown. The explicit calculation of statistically signifícant correlations from raw 
data requires the processing of large dataseis and Fourier or power spectrum analysis to extract the 
dominant wavelengths. In order to avoid this computational burden and improve the flexibility of the 
analysis scheme, the OMEGA OI scheme represents the covariances by an ideal ized autocorrelation 
function. One common correlation function, used by DiMego (1988) among others is: 


A* = 


O* exp- {~r ^ 

_L L 


1 + k p ln' 


V 

\ Pk J 


(189) 


where D[k = distance between points 1 and k; 

L = a characteristic length scale; 

kp = constant which determines the sharpness of the vertical correlation function; 
pi and pk = pressures at points 1 and k. 


In a 2-dimensional OI scheme, kp would be set to zero. The factor C] W serves to reduce the correlation 
between points when one is located over land and the other is located over water. C] w is set to one for 
grid cells that are not within 100 mb of the surface, when the two points are both located over the same 
surface type, and for the SST analysis. Otherwise it is set to 0.33 for grid cells at the surface an 
increases linearly to 1.0 for cells 100 mb above the surface. 

Since most operational forecast systems derive their first-guess from the same source and are integrated 
on the same grid each day, some reasonable estímate can be made of £, the ratio of the observational 
error to the first-guess error. OMEGA derives its first-guess from any number of sources and is run on 
grids located anywhere on the globe with grid point distributions that can vary from run to run. For this 
reason, the estimation of £ is considerably more difficult. Currently £ is assumed to be 0.3 (0.37) for 
rawinsonde (surface) height errors, 0.45 (0.9) for rawinsonde (surface) wind errors, 0.45 (0.75) for 
rawinsonde (surface) relative humidity errors, 0.3 (0.75) for rawinsonde (surface) temperature errors, 
and 0.3 for SST errors. 


The multivariate algorithm in the OMEGA data preprocessor is based on that described by DiMego 
(1988). The multivariate problem is similar to the univariate problem with the exception that correlation 
functions between different variables (geopotential, u and v) are multiplied by an additional factor. The 
factor is based on the assumption that the geostrophic relation applies in the coupling of the wind and 
mass fields. Of course, in the tropics, within the PBL, in the presence of strongly curved and/or 
unbalanced flow, this may not be entirely true. But the approximation works fairly well in most cases 
with the addition of a coupling factor G that reduces the coupling near the surface, near the top of the 
model domain and at low latitudes. The user has the option of performing a univariate analysis if a high 
resolution analysis in regions of highly unbalanced flow. The additional factors as described in DiMego 
(1988) are: 
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where: 

subscripts u, v and z represent u wind component, v wind component and geopotential 
observations; 

puz» pvz» p Ul u v pv,v 2 , and p uv represent the correlation factors between variables of different types; 
lat and Ion = latitudes and longitudes (in radians) at the location of the observation or gridded 
valué represented by the subscript; 
i = the analysis grid point location; 

G = degree of geostrophic coupling between the mass and wind fields; and 
r e = radius of the Earth. 


As an example, consider the correlation between a height observation’s deviation from the first-guess 
and a u-wind component observation’s deviation. If the two observations are at the same latitude, one 
would expect no correlation, since a change in the height field at one location does not contain any 
information about a change in the north-south height gradient at that location or any location at the same 
latitude. A positive height perturbation, however, would be correlated with a positive (negative) change 
in the height gradient immediately to its north (south). For this reason, one would expect a positive 
(negative) correlation between height perturbations and u-wind perturbations at locations to the north 
(south) of the height perturbations. 

The geostrophic coupling factor G is set a its máximum user specified valué G m i only at mid and high 
latitudes (more than 25° from the equator) in the mid-troposphere. At mid-latitudes, G increases linearly 
from 0.5G m i at the surface to G m i 175 mb above the surface. Above 250 mb, it begins to decrease 
linearly and reaches a valué of about 0.80G m i at 100 mb. It varíes latitudinally as in Dey and Morone 
(1985). Points greater than 25° latitude from the equator retain their full coupling. The coupling is set 
to zero for points within 10° latitude of the equator. The coupling for points that are between 10° and 
25° from the equator is determined by the following relationship: 

c G.,{cos[M-10)2-180]+i} (195) 

Recommended valúes for G m i are specified in Table 8. 

Table 8. Recommended valúes for user specified OI variables. 


Variable 

Description 

Valué 

Gml 

Máximum geostrophic coupling 

0.85 

L s 

Correlation weighting radius at surface 

450 km 

Lr 

Correlation weighting radius in free atmosphere 

700 km 
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The wind-wind correlations on the diagonal of equation (6.4) are multiplied by an additional factor, H. 

H increases linearly from 0.5 at the surface to 1.0 at 175 mb above the surface. It is then constant up to 
175 mb. Above 175 mb, it begins to increase linearly again and reaches 1.2 at 100 mb. The net effect is 
to decrease the weight of wind observations relative to geopotential observations in the PBL and 
increase their weight relative to geopotential observations in the stratosphere where geopotential errors 
are accumulated from the entire sounding. 


The correlation weighting radius, the length scale in Equation (189), is 
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where Ls (L r ) is a length scale selected by the user for surface (rawinsonde) observations and Ap is the 
pressure difference between the surface and the grid point. This allows for a smooth transition between 
the relatively high resolution surface data and the coarser rawinsonde data. The weights are entered into 
the system of equations and the matrix is inverted to calcúlate the correction at each grid point. 
Recommended valúes of L$ and L r are listed in Table 8. 


The correlation matrix size is set to 27 at the surface. It decreases approximately linearly to 15 at 
200 mb above the surface. Above this level it is set to 12. The matrix size is always set to a valué that 
is divisible by 3 so that a height, a zonal (E-W) wind and a meridional (N-S) wind observation from 
each location can be used. Twelve spots in the matrix are reserved for rawinsonde observations from the 
4 nearest soundings to the grid point. One u wind component, one v wind component and one 
geopotential observation are included from each sounding. These observations are interpolated to the 
pressure level of the analysis grid point, unless the analysis pressure is within a data void región for a 
given sounding. In this case, the nearest observation in the sounding is chosen. For the purposes of 
Equation (189), the observation pressure is assumed to be that of the nearest actual observation in the 
sounding. The remaining spots in the matrix are reserved for observations from the nearest surface 
stations. Surface pressure observation increment hydrostatically converted to a height increment and 
wind competence are selected from each station. At the surface, this leaves room for the five nearest 
surface stations, while there is only room for one 200 mb above the surface and none above that. In the 
univariate case, the matrix is filled with observations of the same type, so that 3 times as many 
observing stations can be used. 

The analysis proceeds horizontally, layer by layer, cell by cell. The analysis is performed at the top 
center of each OMEGA grid cell as well as at the bottom of the lowest layer of cells. Corrections to the 
first-guess field for height and wind are computed using the multivariate algorithm while those for 
relative humidity are computed using the simpler univariate algorithm. If the user sets G m i to zero, all 
corrections are computed with the univariate algorithm. Corrections are computed for an entire layer 
and smoothed with a 4 point smoother before being added to the first-guess. The smoother sets the 
valué at a given cell to 50% of the original valué plus 1/6 the valué at each of the 3 adjacent cells. The 
height corrections are converted to pressure corrections since OMEGA defines grid locations at constant 
heights. The variables are then vertically interpolated to cell centroids. Temperature is analyzed only at 
the surface. At cell centroids, it is derived hydrostatically from the pressure at the top and bottom of the 
cell. If a level is less than 10 mb from the last analyzed level, the analysis at this level is skipped. The 
corrections to the first-guess are then interpolated from the corrections at the nearest analyzed layers 
above and below. This reduces the problem of noisy temperature fields due to small differences in the 
height correction at closely spaced analysis layers. 
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The surface temperature data is used to hydrostatically adjust the surface height increment in the lowest 
levels as in DiMego (1988). At the ground, the surface height increment will be used directly. In the 
first few model levels above the ground (up to about 40 mb above the surface), the surface height 
increment is hydrostatically corrected by assuming that the surface temperature perturbation is valid 
over this layer. The corrected height increment is considered to be valid at the model level. In effect, a 
mini-sounding is created from each surface observation. The profile extends about 40 mb above the 
ground. Only the nearest height increment to a given level is included in the analysis matrix at any 
given grid point. This has the effect of influencing the thickness field near the ground so that the 
hydrostatically derived temperature field closely resembles a direct surface temperature analysis. In 
summary, the procedure to obtain an estimate of the analysis variable Aj with the OI scheme is: 

(1) determine the matrix size and the observations which go into it as described above; 

(2) calcúlate the autocorrelation function valué for each pair of observations and for each of grid 
point location-observation pair using Equations 189 through 195; 

(3) solve equation set 6.4 using a standard matrix inversión technique (Gaussian elimination), to 
obtain a weighting coefficient (wk) for each observation; and 

(4) substitute the wk valúes into Equation 186 to obtain an estimate of Ai at the grid point. 

6.6.5 Adjustment of Analyzed Valúes. 

The final process in the preparation of the initialization file is the adjustment of the vertical temperature 
profile. Since height rather than temperatures are analyzed by the 3-D OI scheme, temperature profiles 
will sometimes inelude regions of superadiabatic lapse rates or temperature spikes. This is particularly 
trae near the surface where closely spaced analysis layers mean that a small error in the height analysis 
at one level can produce large temperature anomalies in the layers above and below. Spikes are 
removed, one at a time, starting with the most significant spike. The spike is removed by adjusting the 
height of the layer interface which lies between the grid point that contains the spike and the adj acent 
grid point which most resembles a spike of the opposite sign. This method is based on the assumption 
that a small error in the height analysis at one layer interface relative to the interfaces above and below 
will introduce a warm-cold dipole in the temperature field at the layers immediately above and below. 
Once all significant spikes are removed, superadiabatic layers are removed in a similar manner. Often, 
the removal of temperature spikes also removes superadiabatic layers so that this final step is not 
needed. The adjustment of the intermediate height is done so that temperature lapse rate between the 
grid point above and below is set to the mean lapse rate of a layer which extends one level above and 
below any contiguous layers which contain a spike or superadiabatic lapse rate. 

The lapse rate adjustment is an iterative process in which the largest temperature spikes and the most 
superadiabatic lapse rates are removed first. The process continúes until no superadiabatic layers remain 
and the largest temperature spike is below a threshold magnitude. Layers near the surface with a 
modérate superadiabatic lapse rate are allowed to remain. The interpolation of variables from cell tops 
and bottoms to centroids as well as the hydrostatic calculation of temperatures from the pressure field 
and the interpolation of temperature to the new pressures at the cell centroids which result from the 
surface pressure analysis are all an integral part of the lapse rate adjustment process. 

6.6.6 Generation of Analysis Performance Statistics. 

At this point, temperature, humidity, winds and heights are analyzed to the mandatory pressure levels at 
rawinsonde sites and to the locations of surface observations for verification purposes. A summary of 
the RMS differences between the 3-D dataset both before and after the update by the 3-D OI scheme are 
written to the "ANAinfo" file. 
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6.6.7 Final Processing and Writing the "INIT" File. 

Next, the valúes in all of the arrays that contain variables which define the properties of the earth’s 
surface are checked. If any of these arrays contain "flag" valúes because no actual data was available, a 
default valué is assigned to that variable. The default valúes for each variable are deterrruned by the 
valúes specified by the user in the "prepomega.opt" file. 

The output from subroutine ANALYSIS consists of an "INIT" file in which the 3-dimensional variables 
and the skin temperature have been updated with the observations. 

6.7 Preparation of Lateral Boundary Conditions. 

The lateral boundary condition data for the OMEGA model is orchestrated by the RUNPREP script. It 
takes advantage of the capabilities of PREPDAT and GENINIT to be run in both the Initialization mode 
described above and Boundary condition mode. A sepárate BC file is created for each boundary 
condition time. The initial file "BC.tagO" is created when PREPDAT and GENINIT are run in 
Initialization mode. RUNPREP then runs the two modules once for each subsequent boundary 
condition time to create the files "BC.tagl", "BC.tag2" and so on. 

6.7.1 Execution of PREPDAT in Boundary Condition Mode. 

When PREPDAT is run in Boundary condition mode, no observational or surface property data is 
processed. A gridded data file is read in from one of the gridded data sources listed in Table 7. The data 
is processed in the usual way, but it is interpolated horizontally only to the boundary points that are 
created by reflecting each cell that borders on the edge of the domain across the domain boundary. The 
processed data is then written to the "GUESS" file. 

6.7.2 The Execution of GENINIT in Boundary Condition Mode. 

When GENINIT is run in Boundary condition mode, it reads in the "GUESS" file written by 
PREPDAT’s Boundary condition mode execution. It then interpolates this data vertically from pressure 
surfaces to the cell centroids of the boundary cells. It then writes out a "BC" file identical in format to 
the "BC.tagO" file without writing an "INIT" file. RUNPREP renames this file to "BC.tag#" where 
number represents the sequence of the file in the list of boundary condition times. 
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Section 7 


Model Boundary Conditions 


7.1 Lower Boundary Conditions. 

The OMEGA model prognostic variables ( pu , pv, pw, p0, andpq ) at the lower boundary are treated in 
the following form 

A pC = P av M + ta V e*P> (197) 

where C, represents u, v, w, 0, and q. 


The lower boundary conditions for and p (represented as g) are formulated with the aid of Monin- 
Obukhov similarity theory for the atmospheric surface layer. This theory allows one to derive 
relationships between the valúes of ground level prognostic variables at the first ( g s ), and second 

computational level ( g 2 ), and a diagnostic level (£,) 

f, =í ! í > .+(l-/ > .k„ (198) 

where 
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(199) 


and y/ m and y/ h are universal functions of nondimensional height £ = zlL. The revised formulas for 
\¡/ m and if/ h proposed by Beljaars and Holtslag (1991) are used. The two lowest computational levels of 
the model are assumed to be placed within the atmospheric surface layer. 


At the ground surface, a no-slip condition is imposed for the wind so that 

m(z 0 ) = v(z 0 )= w(z 0 ) = 0. (200) 

The hydrostatic equation is used to obtain the surface pressure at level z 0 (the sub-terrain “ghost” cell 
layer) from level 2 (the lowest atmospheric layer) OMEGA pressure, while a zero gradient boundary 
condition is imposed on density, p. For all Eulerian tracers (Q t ) zero flux boundary condition (i.e., 
complete reflection) 

K^%- = 0 ( 201 ) 

ai 

are assumed at the ground, as the surface absorption rates are not well known. Lower boundary 
conditions for turbulent kinetic energy, e, is taken to be 

de 

^• = 0 . ( 202 ) 

oz 


87 



7.2 Lateral Boundary Conditions. 

While OMEGA simulates processes that occur over all of the atmospheric System, usually only a portion 
of the atmosphere is chosen for simulation. This requires certain assumptions to be made regarding the 
portion of the physical domain that is outside the computational domain. These assumptions are applied 
as boundary conditions to the edge of the computational domain. In OMEGA, the bottom boundary is 
the surface of the earth, either land (with varying surface properties) or water. The lateral and top 
boundaries are continuations of the atmosphere itself. Each of these boundaries have to be treated in an 
appropriate manner. 

The lateral boundaries are open and should allow the unimpeded flow of air (and associated properties) 
through them. Outflow boundaries are deterministic, while the inflow boundaries are not. In cloud scale 
models, in which the time scale is small enough to assume that the environment outside the model 
domain is unchanged during the simulation, the inflow boundaries are usually prescribed from 
environmental conditions. In mesoscale and larger models this assumption can no longer be made. The 
environment outside the computational domain has to be derived from larger scale forecast tools such as 
the NCEP ETP, Nested Grid Model (NGM), Médium Range Forecast model, (MRF), or the FNMOC 
NOGAPS model. 


The key requirement in connecting the solution within the computational domain with the outside 
environment is that acoustic disturbances be allowed to propágate out with minimal reflection. To 
ensure this we use a radiative boundary condition based on the method of characteristics (Hedstrom, 
1979). Presently, a 1-D formulation is used, although it can be upgraded to multidimensional treatments 
in the future. Consider the linearized I-D gas dynamics equation in matrix form: 
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and a local coordínate system is used such that the boundary flow, Vb, is positive when directed 
outward. There are three characteristic speeds, X¡, derived from the eigenvalues of A: 


(204) 


X 0 = Vb, X+ = Vb+c, and X. = Vb-c. 


(205) 


where c is the local sound speed. The characteristic equations along dx/dt=A,¡ are 
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where 1¡ are eigenvectors of A: 
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From (204), we get: 


dp=dp/c 2 along dx/dt = Vb, 


(208) 


±pc dVb+dp=0 along dx/dt = Vb±c, (209) 

To minimize reflection, signáis along the incoming characteristics should be “quiescent,” i.e., 

1, • dU/dt = 0 for X. = Vb-c, and for X 0 = Vb when Vb<0 (inflow). We can now write the boundary 
condition in the form of initial valué equations for the boundary variables Vb, p, and p. For an inflow 
boundary, Vb<0, we get 
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where —— indicates that the gradient should be evaluated inside the computational domain. For the 


outflow boundary, we get 
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Ap = —r A P + ^~ Uto Uto 


The three boundary variables are advanced in time explicitly. Finally, at the lateral boundaries, zero 
gradient boundary condition is applied for turbulent kinetic energy. 


7.3 Upper Boundary Conditions. 


Since the top boundary is also assumed to be open, an appropriate radiative boundary condition should 
be imposed on it, which will allow vertically propagating gravity waves to escape the model domain 
without being reflected back into it. This is more complicated than for sound waves. Klemp and Durran 
(1983) show that when the gravity wave phase speed is larger than the mean horizontal wind, the vertical 
phase speed and group velocity are in opposite directions, i.e., a simple characteristic treatment will not 
work. They formulated a top boundary condition scheme that can radíate gravity waves but is both 
difficult and computationally expensive to implement on an unstructured grid. Henee in OMEGA we 
are using a Perkey-Kreitzberg (1976) type sponge defined over the top few layers. This is in addition to 
the acoustic radiative boundary condition described in the previous section. Since the phase speeds for 
sound and for the gravity waves differ by an order of magnitude, the acoustic boundary condition should 
have no influence on gravity waves. 
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The sponge treatment should be adequate in most cases especially if the sponge layers are placed well 
above the tropopause. According to this scheme the prediction of a dependent variable 0can be written 
as 


0 'o +A '(r)= 0 'o(r) + 
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where subscripts m and Is refer to the model predicted tendency and a large scale tendency which is 
derived from a larger scale forecast model (LFM or NGM models), and w(r) is a weighting function 
defined as follows: 



for r > r. 


sp 


for r<r. 


sp 


(213) 


where r sp specifies the start of the sponge layer, and a is usually set at 0.5. Thus the valué at the 
boundary is completely prescribed by the large-scale flow and the valúes at the bottom of the sponge 
layer are the model calculated valúes. If the large-scale tendency is not used this scheme becomes a 
porous sponge in which the boundary is not tied down to any specific valué but allowed to float in 
response to the interior flow. 


Finally, turbulent kinetic energy, e, vanishes at the top boundary (i.e., e=0). 
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Section 8 


Atmospheric Dispersión Model Structure 


8.1 OverView. 


The Atmospheric Dispersión Model (ADM) is embedded within OMEGA. The model is composed of 
Eulerian and Lagrangian transport algorithms that advect massless tracers, gases, and/or solid and liquid 
particles using the OMEGA-resolved wind field while simultaneously solving a diffusion model that 
simulates the effect of unresolved physics including turbulence. Massless tracers and gases are in 
thermodynamic and dynamic equilibrium with the ambient atmosphere whereas solid and liquid particles 
are allowed to settle out with a size and density-dependent vertical slip velocity due to gravity. The 
Lagrangian routines provide a means of modeling point releases. The Eulerian tracer capability provides 
a means to track releases of gaseous species originating from OMEGA grid cells where it is assumed 
that the released gas instantaneously and uniformly filis the OMEGA cell containing the source. The 
Eulerian tracer capability also provides a first order comparison against the Lagrangian particle transport. 

8.2 Eulerian Model. 


An Eulerian tracer methodology is applicable to large area releases. It is assumed that the released gas 
instantaneously and uniformly filis the OMEGA cell containing the source and that the source strength 
(release rate), remains constant through the OMEGA run. The OMEGA Eulerian tracers are currently 
assumed to be massless and henee are advected and diffused with the ambient air. The conservation 
equation for the tracer field can be written as 


dQ t 

dt 


= -V(0V) + S, 


(214) 


where Q, is the tracer concentration, S is the subgrid, and S, represents the tracer sources/sinks including 
subgrid transport fluxes. 


Sources can be specified via their source strength, latitude, longitude, and altitude locations. The source 
specification can be allocated through the OMEGA graphical user interface or by editing the OMEGA 
configuration file. An example of an Eulerian release is shown in Figure 27. 


8.3 Lagrangian Particle Model. 


The ADM Lagrangian particle transport provides a comprehensive injection and dispersión capability for 
puffs or discrete mass elements. The user can choose to initialize múltiple release locations for either 
massless or settling particles. In addition, these particles can represent individual masses or the 
centroids of diffusing puffs. The user also has the ability to predetermine the altitude of release, the 
injection time interval, the start and stop time of the injection, the number of particles (or puff centroids) 
initialized per injection time interval, and a number of particle physical parameters. The resulting 
particle output can be post-processed to produce concentration fields for comparison to observed data or 
simply as a means to visualize OMEGA flowfields. 
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8.3.1 Particle Transport. 

The ADM Lagrangian model uses a particle-in-cell technique to simúlate the transport of particles in the 
OMEGA wind field. Particles are advected using the OMEGA-wind components. In addition, particle 
velocity is adjusted so that particles are allowed to fall out with a vertical slip velocity due to gravity. 


In order to transport a particle horizontally and vertically in response to the OMEGA-generated wind 
field, ADM first averages the domain cell-centered spherical (u, v, w) components of velocity to the 
domain cell vértices. Next, the velocity components are interpolated to the projection of the particle 
position in the model layers both above and below the actual vertical location of the particle. Finally, a 
bilinear interpolation is performed between the model layers above and below the vertical location of the 
particle to get the components of velocity at the particle location. 


Massless particles and the centroids of vapor puffs are advected with the components of velocity 
generated in this manner. Since they are in thermodynamic and dynamic equilibrium with the ambient 
atmosphere, they follow an average of the local wind flow exactly. Massive particles on the other hand, 
must be adjusted for particle slip velocity due to gravity. The slip velocity is determined by calculating 
the terminal velocity, V¡, using modified Stokes drag. The terminal velocity calculation starts with the 
particle momentum equation 
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Setting 


w = the component of particle velocity perpendicular to the ground, 
a = drag coefficient multiplier for rough spheres, 
g = acceleration due to gravity, 

R e = Reynold’s number, 

A p = cross-sectional area of the spherical particle, and 

p a = ambient air density. 

dw . 
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Particles Crossing the top and side boundaries are removed and the resulting particle arrays compressed. 
Those particles impacting the surface can either deposit or perfectly reflect. 


8.3.2 Particle Tracing Methodology. 

Once particles or puff centroids are initialized, their origin and subsequent locations on the OMEGA 
unstructured grid are updated in the vertical using a simple comparative search on particle and layer 
altitudes. Horizontal locations are updated using a successive nearest-neighbor search routine based on 
the use of vector cross producís (Lottati, 1993). 
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The basic steps of this routine are illustrated in Figures 27 and 28. In Figure 27, a partióle, AH, the 
relative horizontal location of one element of the unstructured grid, is evaluated with respect to each 
vertex location of the element. The search routine computes the cross product between the vector 
formed from the first vertex location to the current particle location, A, and the vector formed from the 
first vertex location to the second vertex location, B. A positive valué of AxB (positive area with 
respect to B), indicates that the particle is located in the element with respect to the first vertex location. 
Additional cross producís are computed with respect to the other two vértices with the result that three 
positive retums indicate that the particle is located in the element. 

In Figure 28, the particle, now at location X2, has advected some distance from its original location to a 
point out of the original element. The search routine again evaluates the cross producís using the oíd 

vertex locations and the new particle location. Evaluation of the cross product retums a negative AxB 
(negative area with respect to B ), for the first vertex indicating that the particle has crossed the element 
boundary and now possibly resides in the adjacent element. Stored grid information updates the element 
identifier on the right hand side of the original element boundary, and the search routine moves on to the 
adjacent element for testing. If the particle had advected several elements away, the scheme would 
continué to advance through neighboring elements in a directional sequence dictated by which vector 
cross producís retum negative valúes. The search would continué until three positive cross producís are 
retumed for one element, thus indicating that the particle has been located and its new cell id 
determined. 

8.3.3 Particle Dispersión. 

ADM has the capability to simúlate the dispersión of non-buoyant particles using either a Lagrangian 
particle model or a formulation for puff diffusion. The particle model treats each particle as a discrete 
element of mass. This differs from the puff model which treats each particle as the centroid of an 
expanding puff. The Lagrangian particle model simulates the dispersión of pollutants or by tracking a 
large set of Lagrangian particles moving, at each time step, with pseudo-velocities. These pseudo- 
velocities simúlate the effects of the two basic dispersión components: (1) transport due to the mean 
wind, and (2) diffusion due to turbulent velocity fluctuations. Subsequent positions of each particle (x, 
y, z) representing a discrete element of pollutant mass are computed from the following 


x(t + At) = x(t)+[u(t) + u'(t)]At , 

(217) 

y(t + At) = y(t) + [v(í) + v\t)]At , 

(218) 


and z{t + At) = z(í)+[w(t)+w'(r)]Ar (219) 

where ü , v and vv are the mean wind components obtained directly from OMEGA model, and u’, v’ and 
w’ are the corresponding subgrid scale turbulent velocity fluctuations whose statistics are derived from 
the boundary layer formulations used in the OMEGA model. 

The particle turbulent velocity fluctuations are calculated using a Markov chain scheme 

u'(t + A t) = R (At) u'(t) + u"(t + At) 


v'(t + A t) = R ,(At) v'(t) + v"(t + At) 


( 220 ) 

( 221 ) 
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Figure 27. A partióle, XI, located in one element of the OMEGA unstructured grid, is 

evaluated using AxB. The cross product between the vector formed from the first 
vertex location to the current partióle location, A , and the vector formed from the 
first vertex location to the second vertex location, B indicates a positive valué of 
area with respect to B. 


A 



Figure 28. A partióle, XI, located in one element of the OMEGA unstructured grid, is 

evaluated using AxB. The cross product between the vector formed from the first 
vertex location to the current partióle location, A , and the vector formed from the 
first vertex location to the second vertex location, B indicates a negative valué of 
area with respect to B. 


w'(t + A t) = /? W .(At) w'(r) + w"(t + A t) 


( 222 ) 
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where u” v” and w” are purely Tandom uncorrelated turbulent velocity components. Distributions of 
these random turbulent velocity components are Gaussian with zero mean valúes, i.e., they are 
completely characterized by the following variances 


<=<(>-<) 

a 2 = ct 2 (l -R 2 ) 

r w . w'/ 


( 223 ) 


(224) 


(225) 


where R u ,, R^,- and R w . are auto-correlation coefficients for each velocity component and assumed to be 
exponential 
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(228) 


where T“ , T[ and T™ are the Lagrangian integral time scales of the sub-grid scale turbulence. 

The Lagrangian turbulent statistics (o u , <y v , o w . Ti, Ti, and Ti ),in each of the three component 
directions must be determined. An important advantage of Turbulent Kinetic Energy (TKE) closure is 
that it allows one to determine a variety of turbulent variables including the Lagrangian turbulent 
statistics (Uliasz, 1990). The variances of wind velocity components used in the Lagrangian particle 
model can be written in the following form 


^.=2el~2A t (S„G„ +S H G„) + 6A,S„ j 

íi f 2 s Y 

al=2e\--2A,{S„G„ +S„G H ) + 6A l S„ — ■ 

<7 W = 2ej— — 2A, ( S m G m + S h G h ) + 6 A¡S H G H j 


(229) 


(230) 


(231) 


In addition, the Lagrangian integral time scales are related to the diagonal elements of the diffusivity 
tensor, Kxx, Ryy, Rzz, 
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where K xx =í4l¡S xx , K yy = i42~eS yy , = t4leS z 

with the functions Sjj expressed by the following parameter. 
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8.4 Probabilistic Puff Model. 

Continuous plumes are commonly represented by subdividing them into discrete, overlapping volumes, 
often called “puffs”. The major components of the probabilistic puff model are discussed in detail 
below. 

8.4.1 Generation of Puffs. 

The model generates a new group of Lagrangian puffs (N) at a user’s defined timé interval (A/). The 
number of puffs released, N, depends on local wind speed and turbulent diffusion at the source. N is 
chosen in a such way that the distance between 2 puffs should be approximately equal to cr h / 2. Each 
puff generated represents a discrete element of the total material released into atmosphere. The amount 
of the material, Q, is the same in each puff so that 

0 = 9— (237) 

N 

where q is the source strength (kg/s). 

8.4.2 Motion and Modification of Puffs. 

All motions and changes in size are determined by conditions at the puff center by the following three 
operations: (a) the transport of puff centroid by the mean wind, (b) the distribution of material within 
each puff (around the puff centroid) in space according to a Gaussian function whose <7 X , < 7 V , and <r z 
valúes increase with travel time or distance, and (c) plume rise. 

Large scale air modification of each puff is discussed in detail in Section 8.3.1 (particle transport by the 
mean wind). For puff diffusion, the standard deviations of the Gaussian distribution are estimated based 
on Taylor’s homogeneous diffusion theory which use turbulence characteristics from the turbulence 
closure scheme used in PBL model 

cr x (t + At) = G x (t) + <J a At for/<2r L “ 
and (238) 

a ](if + Ai) = a x (/) + 2T u L a]At for t > 2T L 
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where the variance of wind velocity component, o u , and the Lagrangian integral time scale, T[, are 
discussed in Section 8.3 (Lagrangian Particle Model). In a similar way a and a z are determined 
(Uliasz, 1990). 


For the plume rise treatment, the work of Briggs (1969,1979) provides the basis. The following plume 
rise equations are applied separately to each puff centroid, using the wind and stability valúes at the 
point. This allows spatial and temporal variations of atmospheric stability to be treated. For plume rise, 
the goveming parameter is the buoyancy flux F(mV 3 ), given by 


gv,d)(T,-T) 


(239) 


where g gravitational acceleration, v e stack gas exit speed, d t stack exit diameter, T e stack gas exit 
temperature, and T ambient temperature. Then, four different equations (depending on wind speed and 
stability) are used to determine plume rise. For purpose of determining plume rise, stability is needed 
and defined by 


S = !L— 


(240) 


The equations for plume rise w when S > 4.2 x 10 VIH (stable conditions ) are 
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where H s is stack height, V is horizontal wind speed, and V s is critical wind speed separating ‘calm’ 
from ‘not calm’ conditions for stable conditions and is defined by 


V s = 0.14F' S 


(243) 


Note that at low wind speeds, travel time is used to determine plume rise instead of travel distance. 
Máximum plume rise under stable conditions is achieved when 

t = l0 TT (244) 

V S 

A similar approach is used in the following equations when S > 4.2 x 10 -2 V 2 / H 2 (neutral and unstable 
conditions) 




when V >V_ 


(245) 




when V <V. 


(246) 
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where V n is critical wind speed separating ‘calm’ from ‘not calm’ conditions for neutral or unstable 
conditions and is defined by 

v. =-!^—(lOtf.f 3 (247) 

where is the máximum plume rise allowed during light winds and neutral or unstable conditions. It 
is set to 3000 m to ensure that all parameters remain fínite, but that the plume does not impinge at 
ground level. 

8.4.3 Puff Splitting Scheme. 

As the puffs grows, the local representation of the turbulence and velocity fields using puff centroid 
location becomes increasingly inaccurate. When the meteorological fields are inhomogeneous, the 
accuracy of the calculation can only be maintained by splitting puffs into smaller components which can 
sample the variations in the meteorology explicitly. Puff splitting scheme is implemented in such a way 
that when the growth of a puff extends either vertically or horizontally over one cell area of OMEGA 
model, then the puff is divided into two or more puffs in which the sum of their masses equals the mass 
of the initial puff. 

The horizontal growth of a puff continúes until ratio of the puff radius, R, to the grid spacing, G, is 
greater than 1. At this point the puff is divided into four new puffs, each with 0.25 of the initial mass. 
The new positions are orthogonal and equidistant from the oíd position by 0.5 R, with a <7 h valué of 0.5 
of the initial valué but with the same o í valúes. 

In the vertical dimensión this process can be described as follows. When the growth of the puff is large 
enough to fully encompass at least two vertical layers, then the puff is divided into two puffs of mass. 
The new puff centroids have the same horizontal position but have an elevation equal to the layer 
midpoint. The a h valué remains the same but the new valué is proportional exactly the same way as 
the mass. 

8.4.4 Puff Merging and Purging Scheme. 

The key element to the practical application is the elimination of duplícate puffs. The generation of new 
puffs, especially due to vertical diffusion can quickly overwhelm even the fastest and largest computers. 

Therefore, overlapping puffs with Gaussian concentration distributions is merged when they are 
separated by about <7 h or less. This reduces the number of puffs used in probabilistic puff model. Every 
ten minutes all adjacent puffs are checked to see if puff centroids are within o h of each other. If so, the 
material from both is merged into a single puff. The centroid of the merged puff is assigned new 
coordinates and a h valúes that are averages of those in the original puffs. 

Puffs are also removed from further considerations when the puff centroid point is outside of the 
modeling domain. This approach should cause no problems within the area of interest, except when a 
wind reversal brings puffs back into the modeling domain. 
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8.4.5 Concentration Calculation. 

The concentration within each puff is distributed in space according to a Gaussian function whose a x , 
a y , and o z valúes increase with travel time or distance. Knowledge of puff location, the material in it, 
and its geometrical dimensions allows concentration to be determined anywhere by summing 
contributions from each puff. Local concentration valúes are receptor oriented and based on summing 
contributions from each puff using generalized Gaussian distribution with reflection terms from the 
ground and the mixing layer height 


C(x,y,z,t) =■ 
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where C is the concentration, Q is the mass of puff, x,y,z are the position of receptor, x c ,y c ,z c are the 
position of the puff center point, f g is the fraction of cloud mass reaching the ground reflected back into 
cloud mass, f m is the fraction of cloud mass reaching the mixing layer reflected back into the cloud, z m 
is the mixed layer height. 


8.4.6 Surface Dosage. 

The integrated surface dosage can then be defined as 

í 

X(x,y,t) = jc(x,y,Q,t')dt’ (249) 

o 

This quantity is accumulated on a stationary receptor network. 


8.4.7 Calculation of Concentration Fluctuation. 

Time series of all atmospheric boundary layer variables exhibit turbulent fluctuations and air pollutant 
concentrations are no exception. Observations have shown that concentration fluctuations are 
sometimes at least as large as the mean. Therefore, there is a growing need for methods of predicting 
fluctuations of concentration. For example, accidentally released some gases are toxic and their 
flammability is determined by short term concentration levels (e.g., the accidental release of MIC gas 
during Bhopal gas accident). In this emergency case, the concentration fluctuations are crucial to 
predict. On the other hand, the magnitude of turbulent fluctuations of concentration determines the 
uncertainty in air quality models. Uncertainty calculation is also important in most emergency hazard 
predictions. 


Wilson et al. (1982a, b) were concemed with simulating concentration fluctuations from a continuous 
source observed in a wind tunnel and proposed a simplified empirical model based on the standard 
Gaussian formula. This model is used in probabilistic puff model for concentration fluctuations 



( 250 ) 
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where q is the effective variance source strength, b is a constant determined by the form of the vertical 
distribution, U is the effective transport speed, d y and d l are lateral and vertical plume half-widths, and 
G and F are dimensionless functions 
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where the source location parameter, /?, and the sink strength, a, are both assumed to equal 0.6 based on 
analyses of the wind tunnel data, and the constant A = 0.693. 


8.5 Lagrangian Model Modes of Operation. 

8.5.1 Tracer Mode. 

Operating in tracer mode, the ADM Lagrangian model allows the user to inject non-buoyant, passive 
particles into the evolving OMEGA wind field. The user has the option to select the number of injection 
points, their locations and altitudes, and the frequency of particle injection. Additional input fields exist 
for particle diameter, density, and mass which allows the user to physically-characterize massive 
particles. 

Tracer mode particles are assumed to represent discrete mass elements or the centroids of diffusing 
puffs. Particles are removed from the computational domain when they cross the top or lateral 
boundaries. If they impact the surface, massless particles are allowed to perfectly reflect (i.e. no surface 
deposition); massive particle remain on the surface. 

8.5.2 Monitor Mode. 

Monitor mode is designed to provide a monitoring capability over time at múltiple locations within the 
OMEGA domain when an accidental release or spill is expected but the precise release location and/or 
time and/or altitude is not known a priori. The particle output from these monitor station locations can 
then be postprocessed to determine a quick, first-order estímate of the expected atmospheric hazard. 
When combined with air sampling data, monitor mode can also be used to determine the most probable 
geographic source area from which particles have advected. In either application, the user has the option 
to select the number, location, and particle injection interval of the monitor stations. The characteristics 
of the particles (particle type, diameter, and release altitude, etc.) that monitor mode injects into the wind 
field can also be specified. 


100 




Section 9 


Dynamic Grid Adaptation 


9.1 Introduction. 

The accurate solution of any complex computational problem depends on fine spatial discretization of 
the computational domain. The degree of grid refinement is usually constrained by CPU restrictions and 
by run-time expectations. In atmospheric flow simulations, a large amount of CPU time is required by 
planetary boundary layer physics, cloud microphysics and radiation transfer calculations. For real-time 
flow predictions, optimizing the degree of grid refinement, given the CPU constraint becomes crucial. 
Use of unstructured grids simplifies this problem by adapting to regions where a high degree of 
resolution is required, such as shorelines and areas of large terrain gradients (see Section 5). Adapting 
the grid dynamically can make further improvement in solution accuracy and run-time. Solution- 
adaptive grid generation schemes have been widely used in aerospace applications with successful 
results, e.g. Lottati and Eidelman [1994] and Baum and Lohner [1989] and Baum et. al. [1993] have 
applied dynamic grid-generation for simulating flows over aircraft, tracking shock waves, and in 
simulating fluid-structure interactions. Dietachmayer and Droegemeier [1992] and Fiedler et. al. [1998], 
have approached the use of dynamic grid adaptation for atmospheric flows using structured grids. Their 
methodologies have been restricted to idealized problems so far and are not suitable for real atmospheric 
flow simulations, which can inelude complicated terrain features. 

Structured grids are not suitable for dynamic adaptation, because the grid-generation requires a high 
degree of user-interaction and user-expertise. It is unrealistic to have a completely automated grid- 
generation code for real-time predictions. Some atmospheric models such as the Regional Atmospheric 
Modeling System (RAMS) have employed a different approach to overeóme the shortcomings of 
structured grids. In RAMS, for example, there is an option of a finer moving grid within a larger coarser 
domain. One disadvantage of this methodology lies in the fact that the trajectory of the moving grid has 
to be pre-defined and therefore cannot be used for prediction. The other disadvantage of moving grids is 
in errors introduced due to interactions between the coarse and finer nested grids. Using unstructured 
grids eliminates these disadvantages. The main advantage of unstructured grids is the ease with which 
dynamic/solution adaptation can be implemented. There is no longer a need for involved user- 
expertise/interaction for creating topologies of complicated terrain features. The whole procedure can be 
fully automated. Automation of the process is not only highly desirable but can also be required in 
operational settings. The unstructured grid has a smooth and continuous transition from coarse to fine 
regions within the whole domain, thus eliminating interpolation errors, which occur due to interactions 
of different nests in structured grids. 

OMEGA is the only operational atmospheric flow model based on an unstructured grid, and fully 
exploits the advantages and flexibility of unstructured grids. It can adapt dynamically to different 
criteria such as particle location and velocity deformation, etc. This capability is crucial in responding to 
emergeney scenarios such as release of hazardous materials. OMEGA with its dynamic adaptation 
capability has a unique advantage over other atmospheric flow models in providing accurate Solutions 
quickly in an operational setting. 


101 



9.2 Methodology. 


If OMEGA is run with dynamic grid adaptation option, then additional memory needs to be allocated for 
all the arrays. The user needs to define the upper bounds of the size of arrays for memory allocation. 

The user also specifies the adaptation criteria (e.g. particle locations or velocity deformation, etc.) and 
time for invoking adaptation cycle (e.g., every one or two hours). At the time of adaptation OMEGA 
checks if the specified criteria is met. The cells are tagged for refinement or coarsening accordingly. 

The processes of grid refinement and coarsening are discussed in Section 5. Physical variables, which 
are cell-centered, are interpolated to vértices before each adaptation cycle using the pseudo-Laplacian 
weighted averaging scheme. As new cells are created or removed locally a simple linear interpolation is 
done to assign valúes to new vértices and cell centers. All the diagnostic variables are calculated for 
new cells, the particle locations array is updated, underlying terrain is updated, and the mesh is re- 
masked to calcúlate new time steps for the time-splitting scheme. At the beginning of a simulation the 
mínimum elevation is saved. As the grid is refined or coarsened, no cell is allowed to have an elevation 
lower than the initial mínimum elevation valué. This is to insure that no extrapolation is done. The 
original base State variables are also saved and as new cells are created, base State profiles are generated 
for these cells using the initial base State. The whole procedure requires only a few user inputs at the 
beginning of the simulation and after the initial inputs, it is completely automated. 

9.3 Adaptation Criteria. 

Setting the right criteria for adaptation is very important. If the adaptation criteria are particle locations, 
then it becomes a simple step of locating the particles in terms of cell id’s and then tagging those cells 
for refinement. This method can become costly if there are múltiple release sources and the frequency of 
particle release in time is very high. The altemative method is adapting to the flow itself. This is not as 
straightforward as tagging the cells for particle locations. Atmospheric flows are almost always 
turbulent and highly stratified in the vertical. The gradients are weaker (as compared to problems arising 
in aerodynamics, e.g., tracking shock waves) and harder to detect. One advantage in simulating 
atmospheric flows, however, is the large temporal scale of the simulation (usually 12 hours to 48 hours). 
Therefore, even if an event is missed at first, the chances of capturing it at a later time are higher. 
Different adaptation criteria used by OMEGA are discussed below: 

9.3.1 Particles. 

OMEGA was conceived primarily as a hazard prediction model. It already has the unique feature of an 
embedded atmospheric dispersión model. The latest versión of OMEGA also has the capability of 
adapting its grid to particle locations. OMEGA can provide high fidelity Solutions in the regions of a 
particle plume by refining the grid around the particle locations. Coarsening the grid where high 
resolution is not required (e.g., areas from which the plume has already passed), further optimizes the 
problem. The criteria for adaptation is set by defining a Gaussian function around each particle. The 
región cióse to particles therefore gets higher weighting and is tagged for refinement. Areas, which do 
not have particles and consist of relatively fíat terrain are tagged for coarsening. Considerable reduction 
in CPU time is achieved by providing higher resolution only when and where it is needed. 
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OMEGA has the capability to adapt its grid in regions of highly deforming flows. Averaging the valúes 
of deformation in each cell over the first ten levels sets the criteria for adaptation. The máximum valué 
of these averages is calculated which is used to normalize all the average valúes. In this way regions of 
critical flow development are tagged for refinement. 

9.3.3 Adaptation to Other Physical Quantities. 

OMEGA can also adapt its grid to other physical quantities such as absolute vorticity and the potential 
temperature gradient. Work is in progress to define criteria, which is a weighted set of different 
meteorological variables. An important step in this direction would be to create a database of responses 
for different scenarios and use neural networks/artificial intelligence to identify and set adaptation 
criteria. 

9.4 Pseudo-Laplacian Interpolation. 

During the dynamic adaptation, as new cells are generated, all the physical variables are interpolated 
using a pseudo-Laplacian weighted averaging scheme suggested by Holmes and Connell [1989]. 

Rausch[I992] compared the pseudo-Laplacian averaging scheme with area-weighted and reciprocal 
distance schemes and found that the pseudo-Laplacian scheme, although more expensive, 
computationally, yields more accurate results. 

The pseudo-Laplacian for a vertex can be defined as, 
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stability requires the weights, w. to be as cióse to unity as possible. 
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Since, linear functions have zero Laplacians, we can write, 
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Defming the cost function, C as 
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We now have an optimization problem of minimizing the cost function given the constraints in 
Equations (256) and (257). 


Using Lagrange multipliers we can write, 
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9.5 Adaptation to Atmospheric Dispersión. 

Two simulations are compared to illustrate the advantage of dynamic grid adaptation vs. global 
refinement. 

Case A: Global Refinement 

Simulation domain extended from -98.65 degrees to -95.65 degrees in the west-east direction and from 
41.0 degrees to 45.0 degrees in the north-south direction, covering parts of Nebraska, South Dakota, 
Iowa and Minnesota. A variable grid resolution of 5 km to 15 km was specified. Thirty (30) second 
resolution terrain and 5 minute resolution shoreline data was used. The final grid had 2335 cells and 30 
vertical levels, for a total of 70,050 cells (Figure 29, left). In the vertical a stretched grid was used. The 
fírst level near earth’s surface had a thickness of 15 meters and further levels were generated using a 
stretch ratio of 1.15. A particle source was defined in the Southwest comer of the domain with an 
emission rate of 5 particles per minute. The simulation was initialized using the U.S. Navy’s NOGAPS 
gridded data. The simulation was run for a real-time prediction at 15 hours after initialization. 
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Figure 29. Computational grid for Case A simulation (left) and Case B simulation (right). 

Case B: Dynamic Adaptation 

All the parameters in this run were identical to Case A simulation, except the initial grid size. The 
starting grid had 582 cells in the horizontal (Figure 29, right). In this case a variable resolution of 10 km 
to 36 kan was prescribed in the initial grid generated. The adaptation criteria was set to partióle 
locations. 

A refmement cycle was invoked after every hour of simulation time. Cells were tagged by defining a 
Gaussian function around the partióles. Tagged cells were refined and all physical fields were 
interpolated at each adaptation cycle. A coarsening cycle was invoked after every three hours to remove 
triangles from areas where higher resolution was not needed. This was done by identifying the cells, 
which did not contain partióles and were in the regions which were relatively flat and then deleting them 
(Figures 29 and 30). The final grid after 15 hours (Figure 30) into the simulation had 1162 cells. The 
edge lengths varied from 6.3 to 36.7 km. 

Use of dynamic grid-adaptation resulted in more than 40 percent reduction in CPU time and produced 
similar results compared to global refmement. This example clearly demonstrates the inherent 
advantages in using unstructured grids. The simulation required only the initial user inputs (adaptation 
criteria: particle locations, time of refmement: every 1 hour, time of coarsening: every three hours, 
máximum bounds for arrays: 2.5 times the initial grid), after which the whole process was seamless and 
automated. 

9.6 Effects of Improving the Underlying Terrain. 

In each adaptation cycle OMEGA has the option of updating the underlying terrain. The refmement 
cycle resolves finer topographic features, which can drastically improve the solution. In the following 
section three simulations of an idealized case are discussed to examine the effect of growing terrain. By 
introducing an obstacle in uniform cross flow, dynamic pressure builds up around the obstacle and is 
deflected by the obstacle. The perturbation introduced eventually dissipates and the flow reaches steady 
State. Figure 31 shows the time evolution of pressure perturbation field, as the height of the mountain is 
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(a) 5 hours 



(b) 10 hours 



(c) 15 hours. 

Figure 30. Partióle locations for Case A simulation íleft) and Case B simulation (nght). 
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Figure 31. Time evolution of pressure perturbation field for the growing mountain case. Vertical 
cross sections at initialization (left), 8 minutes (center), and 15 minutes (right). 

increased from 350m to lOOOm. Although this is a physical phenomenon, it can create noise in the 
solution. Our results so far have shown that these perturbations do not effect the solution in any adverse 
manner, or they dissipate quickly (on the order of 15-30 minutes). 

Figure 32 (left) shows the pressure perturbation field for the static case. In this case the mountain was at 
a constant height of lOOOm AGL through out the simulation. Figure 32 (right) shows the pressure 
perturbation field for the dynamic case. In this simulation, the mountain was initially 3Ó0m high and was 
raised to lOOOm AGL after +8 minutes into the simulation. 

The effect of improving the resolution of underlying terrain, on the other hand, can be drastic. A simple 
scenario is shown in Figure 33. A single mountain lies in the middle of the domain. All the cells 
surrounding it are fíat. There is a uniform southeasterly flow. Particle locations shown in the figures 
correspond to +3 hours into the simulation. The only difference in the three simulations is the height of 



Figure 32. Pressure perturbation fields for lOOOm static mountain case (left) and 350m- 
lOOOm growing mountain case (right) after 30 minutes into the simulation. 




















the mountain. In Figure 33 (left), the height of the mountain was held constant at 350m AGL. In Figure 
33 (center), the height of the mountain was held constant at lOOOm AGL. In Figure 33 (right), the height 
of the mountain was initially 350m AGL and was raised to lOOOm AGL after +8 minutes into the 
simulation. 

When applied to a full-scale operational run, dynamic adaptation can yield more accurate results by 
providing resolution when and where it is needed. The following figures are for two simulations run for 
the Four Comers area, which is characterized by complex terrain features. The domain and the initial 
conditions for both simulations were identical. The domain extends from -110.55 to -107.61 degrees in 
the west-east direction and 35.61 to 38.61 degrees in the north-south direction. The initial grid had 437 
cells in the first horizontal layer and thirty vertical layers for a total of 13,110 cells. The edge lengths 
vary from 6.81 km to 42.18 km. In the first simulation (Case A), the grid was not refined during the 
simulation. Figure 34 (left) shows the particle locations after +12 hours into the simulation. In the 
second simulation (Case B), the grid was adapted dynamically during the simulation. The adaptation 
criteria was set to particle locations. Figure 34 (right) shows the particle locations for Case B simulation 
after +12 hours. An adaptation cycle was invoked after every two hours for this simulation. The final 
grid had 891 cells in the horizontal layer with 30 vertical layers for a total of 26,730 cells. By resolving 
finer terrain features as the plume evolved, a more accurate representation of the flow field was 
achieved. 

9.7 Adaptation to Velocity Deformation. 

OMEGA is primarily a hazard prediction model and using particle locations as the adaptation criteria 
meets the mission goal. It is desirable, however, to use other physical variables (depending on different 
situations) as adaptation criteria. OMEGA has been configured to adapt to a series of physical variables 
such as velocity deformation and potential temperature gradient, etc. 

A simulation was run with velocity deformation as the adaptation criteria for Southwest U.S. The 
computational domain extended from -112.0 degrees to -96.0 degrees in the west-east direction and 30.0 
degrees to 41.0 degrees in the north-south direction. Terrain data with 30 second resolution and 
land/water data with 5 minutes resolution was used. The initial grid had 2283 cells in the horizontal and 



Figure 33. Particle locations for 350m static mountain simulation (left); lOOOm static simulation 
(center); and 350-1000m growing mountain simulation (right). 
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Figure 34. Simulations for the Four Comers, USA. Static (Case A) simulation results after 12 
hours into the simulation are shown in the leñ panel and dynamically adapting 
simulation (Case B) results after 12 hours into the simulation are shown in the right 
panel. 


This simulation was run for 12 hours. An adaptation cycle was invoked after every two hours. Cells 
with relatively large valúes of total deformation were tagged for refinement. Figure 36 (left) shows the 
final grid at 12 hours into the simulation. The final grid had 3257 cells in the horizontal. As the front 
moved towards the east over Texas and Oklahoma, the grid dynamically adapted to regions of high 
activity. The simulation was able to capture the low pressure system developing over northem Texas 
and Colorado (Figure 36, right). 


Work is in progress to develop new adaptation criteria, as generic weighted functions of different 
physical variables. This way the different complex physical processes occurring in the Earth’s 
atmosphere would be taken into account. Implementation of dynamic adaptation into OMEGA model is 
a major advancement in numerical modeling of the atmosphere. It is the first step, however, and there is 
need for further improvements and optimization of the techniques described in this section. 


thirty vertical levels (Figure 35 left). The vertical grid was stretched. The first vertical level was 15 
meters high and the remaining levels were generated using a stretch ratio of 1.15. The edge lengths 
varied from 20.8 km to 68.9 km. The domain covered parts of Texas, Arizona, New México, Oklahoma, 
Utah, Colorado and Kansas. The simulation was initialized using U.S. Navy’s gridded NOGAPS data 
for 00Z on June 8, 1998. The date was chosen for high frontal activity over parts of Texas, Colorado 
and Kansas. Figure 35 (right) shows the initial sea level pressure in mb. Color scale is from 995.5 mb to 
1017.0 mb, with blue color representing low valúes and red color representing high valúes. 
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Figure 35. Adaptation to velocity deformation. Initial grid (left) and sea level pressure (right) at 
the start of the simulation. 



Figure 36. 


Adaptation to velocity deformation. Final grid (left) and sea level pressure (right) after 
12 hours into the simulation. 

























Section 10 


Parallelization 


OMEGA is CPU bound. Normally less than 1% of the CPU is used for I/O and memory management. 
The computational intensity of OMEGA has been evaluated by the Jumpview profiling system. 
Jumpview determined that OMEGA is dominated by vector operations, which are inherently faster than 
scalar operations; thus, good execution speeds are being achieved. OMEGA exhibits good, but not 
outstanding, vectorization for memory references, primarily due to the indirect addressing inherent in an 
unstructured grid model. OMEGA also showed a computational intensity ratio of flops to memory 
references of 0.99. 

The optimizaron level used in OMEGA is good. It maximizes inlining and vectorization. Inlining is a 
technique to reduce CPU associated with subroutine/function calis. It can cost up to 150 dock cycles to 
cali a single argument subprogram. Vectorization allows a single operation to be performed on "vectors" 
of operands instead of sequentially for each operand. With the optimization level used in OMEGA, the 
overall computational speed achieved is slightly over 40 Mflops. 

Parallelization is a state-of-the-art technique to boost computational speedup by multiprocessing. 
However, unstructured grid models tend to demand significantly more synchronization and geometry 
mapping overheads in theirs parallelization algorithm. Figure 37 shows a speedup projection of 
OMEGA by a power Fortran automatic code-parallelization compiler in terms of factor speedup versus 
processors. The two major reasons for the rather fast diminishing retum: (1) many non-contiguous 
parallelized segments, and (2) large idle time. 

The execution chronology histogram on top of Figure 37 shows that the OMEGA code is translated by 
the power Fortran into many parallelized segments (checkered) and unparallelized segments (striped) 
intermittently interrupting each other. The unparallelized segments represents synchronization bottle 
necks. Another major overhead incurred is inter-processor communication as shown in the short 
checkered bars where the gain in parallelization have been defeated by the expensive cost of message 
passing. 

To overeóme (1) and (2) simultaneously, manual parallelization with a standard message passing 
interface (MPI) is a promising approach. The functionality of MPI’s allows: 

(a) Parallelization in an overarching manner in the OMEGA time advancement loop, making 
message passing minimal 

(b) Pin pointed Solutions to the idle time, or the imbalance load problem, such as techniques for 
accounting the time-masked subeyeling and dynamic domain-decomposition, can be 
incorporated. 

(c) Transportability. 

Figure 38 shows the file structure being implemented. 

Parallelization by domain decomposition of the model domain into contiguous sub-domains is currently 
under development and testing. The decomposition routine is called after dynamic adaptation. Figure 
39 shows an example of the technique applied to the klOm (the Kamisiyah fine mesh about observation 
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Figure 37. Speedup projection of the OMEGA code by a power Fortran compiler. 

points) model domain which is consisted of 4331 cells. It shows how the klOm domain with several 
regions of fine resolution areas is decomposed into 9 and 16 sub-domains resulting in roughly the same 
number of cells in each of the sub-domains. For application on massively parallelized machines the 
scheme is designed to utilize over a hundred processors. 


MPI command cali facilitates the parallelization of the time advancement loop of the model. The test 
case of parallelizing the advancement loop with only the microphysics module tum on has been tested on 
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Figure 38. File structure for the parallelized OMEGA code. 


an IBM-SP2 and an Origin 2000 platforms. Figure 40 shows linear speedups are obtained on both 
platforms for up to 6 processors. 

Sub-linear speedup however is expected for the totality of the parallelized OMEGA code. Interprocessor 
communication is about 10 times more intensive for many of the other modules than that of 
microphysics. For instance, the advection and diffusion modules are among the major interprocessor 
communication overhead consumers due to the evaluations of Laplacian and kinetic energy dissipation 
calculations which require updated neighboring cell information. 

On top of the communication overheads, stacking and unstacking arrays for lumped message passing is 
an inherent overhead for unstructured grid models. These extra CPU costs are both related to message 
passing. However by far the majority of the total overhead is due to uneven load balancing among the 
processors. Figure 41 shows timings of 2 series of runs of a preliminary versión of the parallelized 
OMEGA code for a test case of 1369 cells and 31 variable levels on an IBM SP2 (RS/6000) machine. 
The figure shows timings for the 1-time-mask and the 2-time-masks series. The 1-time-mask series 
means that subcycling of fine mesh cells is not accounted for, henee small time steps of 2.23 s are 
required. The 2-time-mask series does subcycling of the smaller cells and in so doing allows the time 
steps to be relaxed to 4.46 s. The latter is faster than the former for small pools of processors. However, 
for bigger pools of processors, the CPU time gained through the using of larger time step size is 
superseded by the subcycling-induced load imbalance among the processors, i.e., subcycling introduces 
one more level of uneven threadlength across the processors. 
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(a) The klOm model domain topography and 
grid configuration. 


(b) Domain ecomposition of 9 sub-domains 
with roughly equal number of cells in 
each sub-domain. Also visible are the 
boundary cells. 



(c) As in figure (b) except for 16 sub-domains. 


Figure 39. The klOm model with one domain, 9 sub-domains and 16 sub-domains. 


In light of the important contribution of even load balancing to overall parallelization efficiency, there 
exist various techniques to optimize even load balancing. Inverse cell-area weighted domain 
decomposition can be a good way to improve even load balancing, henee, improve the speedup 
performance of the 2-time-mask series of Figure 41. This is an important area of future work of the 
parallelized OMEGA code. For the purpose of illustration. Figure 41 shows the corresponding 
parallelization efficiencies of the runs in Figure 41. 

Another subject of interest is the speedup performance and parallelization efficiency of the OMEGA 
code across different platforms. Workstations are usually more cache limited than mainframe 
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Figure 40. Run time plotted as a function of processor allocation. 

computers. This affects the message passing bandwidth. An investigation of how these are dependent 
on the hardware has been carried out on a cluster of IBM SP2 (RS/6000) machine with 120 MHz clock 
speed and a cluster of SGI Origin 2000 machines with 195 MHz speed. The two platforms are not 
exactly comparable as their Fortran compilers are different. Figure 42, from the same run used in Figure 
41, do show the evolution of speedups and parallelization effíciencies between a mainframe architecture 
and a workstation architecture 



(a) Total Runtime. (b) Parallel Efficiency 


Figure 41. Run time and parallel efficiency plotted as a function of processor allocation. 
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Number o£ Processors Number o£ processors 


(a) Total Runtime. (b) Parallel Efficiency. 

Figure 42. Runtime and parallel efficiency plotted as a function of 

processor allocation for a mainframe and a workstation cluster. 
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Section 11 


XOMEGA Graphical User Interface 


11.1 OverView. 


The OMEGA user interface system is designed around a “case”. An OMEGA case is an encapsulation 
that holds all the information relevant to the simulation. A case ñame is four (4) characters long and 
serves as the ñame of the directory holding all the associated configuration, input, and output files. An 
interactive text-based shell script is used to create and ñame a case. The shell script can also be used to 
examine the status of existing cases. 

The OMEGA System is composed of several distinct modules that are activated in sequence one after 
the another. They are linked together by the interface so that each module gets its specific input and is 
able to direct its output to the next module in line. When creating an OMEGA simulation the user 
makes choices about the geography and meteorology to be simulated. The X Window interface 
(XOMEGA) is the mechanism that translates those choices into input for the other modules. 

After the case has been established, XOMEGA starts. It is used to control parameter settings for the 
other modules of OMEGA. At this stage the cursor and keyboard are used to make selections and input 
data. First, a domain is chosen by using the cursor to define regions on a map. The OMEGA 
computational grid is then generated. The specific meteorology input is chosen for use in creating the 
initial and boundary conditions for the simulation. The meteorological data are fit to the grid by the 
OMEGA preprocessor. Then the simulation is started in a background window. 

11.1.1 X Window Basles. 

XOMEGA is an X Windows Client. An X widows client is an application program that issues requests 
to an X Display Server (user machine) to show text, buttons, and pictures in Windows. These are the 
elements with which the user interaets to make choices about the OMEGA simulation. The mouse is 
used to endose regions on a map, press buttons, and toggle switches; the keyboard is used to input text. 
Several additional Windows will "pop up" and "pop down" during the session. These are called dialog 
boxes and serve to divide the input process according to the modules in OMEGA. Some wamings and 
errors may also appear in dialog boxes. There are a few standard X Windows conventions used in 
XOMEGA listed below. 


These key strokes are used when dialog boxes. The current control is generally outlined in black. 


TAB- 

Shift-TAB- 

Retum- 

OK, Yes, Done- 
Cancel, No - 

Help- 

Clear Form- 
Default- 


Shifts the current control down (or right) one unit. 
Shifts the current control up (or left) one unit. 
Activates the current control, (button press, text input) 
Cióse the current window with affirmative status. 
Cióse the current window with negative status. 

Cancel will clear any text data input in that window. 
Invokes the OMEGA Online Reference System. 
Currently most Help buttons are insensitive. 

Blanks out text input boxes. 

Sets all Controls to the system defaults. 
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Menus are activated by holding down the left mouse button, moving the mouse over the menú, and 
dragging down. The menú Ítems will appear below the title; releasing the mouse button on the desired 
ítem selects that ítem. 

The menus may also be accessed vía keyboard strokes. By holding down the ALT key and pressing the 
underlined letter in the menú ñame the menú will "pulí down". For example, pressing ALT-f the File 
menú will come down. Once a menú is brought down in this manner, the up and down arrow keys are 
used to cycle through the menu’s ítems. Pressing RETURN will active the highlighted item. The left 
and right keys may be used to skip to other menus and "pulí down" on them. Once the menú is down 
the user may release the ALT key. Each of the ítems on the menú also has an underlined letter. 

Pressing the underlined letter will actívate that item. 

At the top of each of the pulí down portions of the menus is a thin dashed line. If the mouse button is 
released on this line the menú will "tear off' from the menú bar and become a sepárate window. This 
allows for semi-permanent access to the menú ítems. The tear off window may be closed by using the 
window manager’s built-in means. Typically this is done by double clicking with the middle mouse 
button on the title of the tear off window. 

Please refer to the operating system manual and window manager manual pages for more information on 
your particular X Window Server. 

11.2 OMEGA Shell Script Interface. 

To invoke OMEGA, the user types "omega" at a UNIX prompt. This starts a shell script that asks a 
series of questions. At any time prior to launching the XOMEGA interface CTRL-C may be used to 
abort the program. There are several sections of the script outlined below. 

11.2.1 Copyright Notice. 

Legal and security notices are given in this section. OMEGA is govemment funded therefore is freely 
distributable to U.S. govemment agencies and their contractors. 

11.2.2 Case Selection. 

The current directory is searched and a listing of all valid OMEGA cases is displayed in tabular form 
(Figure 43). At this stage the user may create a new case by typing in a unique four letter case ID, or 
may view the details of an existing case by typing its case ID. 

The resulting table gives an idea of how much of the case has been worked on. In the left column are 
the ñames of some key files OMEGA uses as the run progresses. The right column shows the status of 
these files, and what they indícate about this case. Figure 44 shows an example of viewing a case. 

The user may continué with the current case by typing its case id. A waming message will alert him that 
some files will be overwritten. The user may view other cases by typing their case ID’s, or the user may 
get a listing of the cases by simply hitting retum. 
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You have selected to view case: zzzz 
Partíal 1isting of /scratch/turner/zzzz : 

case,info Found; The case directory has been created, 

grid.inp Found; The donain has been defined. 

zzzz.adn Missing; The donain has no ALM sources, 

zzzz,grd Found; The grid has been created, 

tine,inp Found; The case has tiñe donain defined, 

net.inp Found; The case has soné net data defined, 

:prep/segl Found; The case has soné net data, 

izzzz.ini Found; The case has an Initialization file, 

¡zzzz,bc00 Found; The case has a BC 0 file, 

Ízzzz,bc01 Found; The case has a BC 1 file, 

zzzz,i Found; The case has a nodel control file, 

zzzzl0000z,out Missing; The case has no analysis output file. 

zzzzlOOOOz.adn Missing; The case has no ADM output file, 

zzzzl0000z.be Missing; The case has no BC output file, 

You nay continué with this case by typing its case id (zzzz), 

Or 

You nay specify a new four character case id to create a new case, 
Or 

You nay type another existing case id to exanine its status, 

Or 

You nay sinply hit <ENTER> to get a 1isting of the OMEGA cases, 
Select case id cr hit <ENTER>: | _ 

Figure 44. View case. 
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11.2.3 Case Creation. 


After choosing a case ñame, a sequence of directories are created that will hold all the files associated 
with the case. First, the case directory itself is created, then two sub-directories are created below that. 
The "grid/" directory will hold input and output files used in generating the case grid. The "prep/" 
directory will hold files used in processing meteorology into initial and boundary conditions. In 
addition, some system configuration files are copied into the case directory, and a few specific 
environment variables are set. All of the creation in this section is automatic. However, there may be 
wamings about trying to re-make existing directories if the user is overwriting a case. 

11.2.4 Launch OMEGA. 

Once the initial case is set up the graphic interface to the OMEGA System begins. XOMEGA launches, 
and now the cursor and keyboard are used to control the system. If the user is not at the machine where 
OMEGA is running, the "DISPLAY" environment variable must be set to the current screen. 

11.3 XOMEGA. 

The XOMEGA application window has four main sections: the Menú Bar, the Operations Column, the 
Status Message Area, and the Drawing Area. Each section has a different purpose and features. The 
Menú Bar spans the top of XOMEGA window and allows for intermedíate and advanced use. The 
Operations Column is the group of buttons at left and is for basic use. Status Messages are displayed in 
a large text box at the bottom. The Drawing Area displays graphics and user selections. Each section of 
the interface is described below. Figure 45 shows XOMEGA as it appears when first opened. 

11.3.1 File Menú. 

The file menú is used to open and cióse files, set preferences, and exit XOMEGA. 

Choosing any of the "Open ..." Ítems pops up a file selection dialog box that allows the user to navigate 
among the directories and files, searching for the desired file to open. The default location is within the 
current case directory. The filter line is set to only show files with the proper ñame and/or extensión for 
the type of file the user is trying to open. 

Once the file opens, the information reads into the appropriate data structure. If the data contains errors, 
an error dialog and the corresponding parameters dialog box pops up. If there are no errors the 
corresponding button in the Operations Column colors gold, to indícate completion of that operation. 

The chart below shows the correspondence between File Menú choice, file type, and Operation button: 

File Menú Choice File Type Operations Button 


Open a Time Input File time.inp 
Open a Grid Input File grid.inp 
Open a ADM Input File adm.inp 
Open a Met Input File met.inp 


Set Time Domain 
Set Spatial Domain 
Set ADM Sources 
Set Meteorology 


Choosing "Preferences" pops up a dialog box (Figure 46) with text entries for the user configurable 
parameters for the Grid Generator and Model. The options and their default valúes are explained below. 
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File Edit Opera tion Window Advanced 


Help 


SOMEGA Opera tions Drawing Area 


Set Time Domain 


mj 


Set Spatial Domain | 
Set ADM Sonrces j 
Rnn GridGenerator ( 


Set Meteojology | 
Feteh Meteorology } 


Ron Pr eproces sor { 
Run OMEGA Model ( 
Exit { 


Operationa! 

Multiscale 

Environment Model with 
Grid 
Adaptivity 


XGMEGA 3.6.0a (C) SAIC 1998-06-21 





XOMEGA Status Messages 


Read in case.info file: 
i Case ñame is aaaa 


Figure 45. XOMEGA. 


Grid Preferences: 

Number of Stretchable Levels (28) - this is the number of levels that are created during the grid 
generation process. 

Vertical Stretching Factor (1.2) - this number determines the spacing of layers in the grid. 

Grid Iterations (10) - this is the máximum number of iterations that the grid generator will 
processes before finishing. 

Please see Section 5 for more information on the OMEGA Grid Generator. 

Model Preferences: 

Field Variable Output Interval (1800 seconds) - the number of seconds between each output file. 
Diagnostic Output Interval (50 steps) - the number of timesteps between each diagnostic report. 

Please see Sections 4,6,7 for more information on the OMEGA Model. 

Choosing "Quit" exits XOMEGA. 
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Grid Preference: 

Number of 
Stretchable Levels 


¡28 


Vertical Stretching f- 


1.2 


Factor 
Grid Iterations íllO 


Model Preference: 


1800 


Field Variable ¡ 
Output Interral (sec) i 

Diagnostic ¡ ^ 
Ouput Interval (steps) j ; 


Clear Form 


Defaults 



Figure 46. Set preferences dialog box. 


11.3.2 Edit Menú. 

The Edit Menú choices offer altemative ways to pop up the various "Set... " dialog boxes. The main 
method uses one of the buttons in the Operations Column. This menú is included for completeness. 

11.3.3 Operation Menú. 

The Operation Menú choices offer altemative ways to execute the extemal programs of OMEGA (Grid 
Generator, Fetch, Preprocessor, and OMEGA Model). The main method uses one of the buttons in the 
Operations Column. This menú is included for completeness. 

11.3.4 Window Menú. 

The only choice on the Window Menú is "XGRED". Choosing this option starts a copy of the XGRID 
program. Section 12 contains more information on XGRID. 

11.3.5 Advanced Menú. 

The only choice on the Advanced Menú is "Set Refinement Disks". This option pops up the Set 
Refinement Disks dialog box. This box is used in combination with the Drawing Area and Map dialog 
to define circular zones within the OMEGA domain. The use of the Refinement Dialog box is discussed 
in Section 11.5.6. 

11.3.6 Help Menú. 

All the choices on this menú actívate the Help dialog box and fill the text window with the appropriate 
help file. The help topics available are: 
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XOMEGA in General 
Setting Time Domain 
Setting Spatial Domain 
Setting ADM sources 
Grid Generation 
Fetching Meteorology 
Preprocessing 
OMEGA Model Execution 

When a topic is selected the appropriate help file is loaded into the dialog’s text window. The text will 
scroll to the top of the section. As the user needs more help sections, they will be appended to the text 
window. The text always scrolls to the top of the most recent section. All Help text remains in the 
window even when it is closed and reopened. The user may refer back to prior help messages at any 
time by simply selecting Open Help from the Help menú. 

11.3.7 Operations Column. 

This group of Controls is presented in order of use (top to bottom). Each button opens either a dialog 
box for detailed input to a module, or a dialog for output from running an externa! module. The buttons 
are listed below with their basic function. The details of each dialog and action, including user 
guidelines for input, are presented in then next few sections (11.4 - 11.11). 


Set Time Domain 

Time Domain dialog box 

Input 

Set Spatial Domain 

Spatial Domain and Map Operations dialogs 

Input 

Set ADM Sources 

ADM Sources and Map Operations dialogs. 

Input 

Run Grid Generator 

Run Grid Generator dialog. 

Output 

Set Meteorology 

Meteorology Parameters dialog. 

Input 

Fetch Meteorology 

Fetch Meteorology dialog. 

Output 

Run Preprocessor 

Run Preprocessor dialog. 

Output 

Run OMEGA 

Run OMEGA dialog. 

Output 


For each button there is an associated dialog box and process. The dialog box pops up when the button 
is pressed. Once the dialog is finished and the process is complete without errors the button tums gold 
colored. If there are errors, or the process is not completed or canceled the button tums gray. 

11.3.8 Status Messages. 

As the user builds the OMEGA simulation messages are received in the text block at the bottom of the 
interface. Errors, wamings, and other status messages are printed in this window. Using the scroll bar 
the user is able to see previous messages. The mouse can be used to increase the size of the Status 
Message window by grabbing the small square on the right hand side above the scroll bar. Be sure to 
drag the window back when using the Drawing Area. 

11.3.9 Drawing Area. 

The Drawing Area is the largest part of the interface and is used to display images, maps, and user input. 
When setting the Spatial Domain, Setting ADM Sources, or Setting Refinement Disks, the mouse is 
positioned in the Drawing Area to "point and click" the input. The behavior of the mouse is controlled 
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by toggle buttons on the Map Operations dialog box. While XOMEGA is drawing, the cursor changes 
from an arrow to a watch. No other activity may take place in the XOMEGA Windows while drawing. 

11.3.10 Map Operations Dialog Box. 

The Map Operations dialog box (Figure 47) is the main control for the function of the cursor when it is 
inside the Drawing Area window. This dialog pops up automatically as needed. There are four drawing 
modes: Zoom Box, Domain Región, Refinement Disk, and ADM Point. The four modes represent the 
geometric shape drawn by the cursor and the input to the interface. 

The mode is changed by clicking one of the diamond shaped buttons on the dialog. The coordinates of 
the cursor are displayed as the user tracks over the map. The cursor will change image for each mode of 
map control. Before drawing in the window make sure that the cursor is the correct image. Use the 
following visual cues to prevent mistakes. 

Mode Cursor Image Drawing Geometry 

Zoom Box four pointed arrow Rectangle 

Domain Región simple cross-hair Rectangle 

Refinement Disk circle Circle 

ADM Point solid dot Point 

The Zoom Box mode is for enlarging sections of the map. The user holds down the left mouse button 
while dragging the cursor into a rectangle shape enclosing the area of the map the user would like to . 
zoom in on. This new portion of the world filis the activity window. 

The Domain Región mode is used to endose speciflc regions of the map that will become part of the 
OMEGA domain. While holding down the left mouse button, drag the cursor into a rectangle shape 
enclosing the area of the map that is the outer-most, or largest, región in the domain. Note that all 
subsequent regions must be within this first región, so make sure it completely endoses the intended 
domain. A total of ten (10) regions may be selected. 

The Refinement Disk mode is used to endose circular areas in the domain that will undergo a 
refinement process in during grid generation. 

The ADM Point mode is used to select points in the domain that will act as sources for the Atmospheric 
Dispersión Model within the OMEGA System. Place the cursor over the desired location and press the 
left mouse button. Up to one-hundred (100) points may be defined. 

There are four buttons on the bottom row of the Map Operations dialog box. The "Cióse" button closes 
the dialog, leaving, the cursor in the current mode, the map in the current State, and with no loss of data. 



Figure 47. Map operations dialog box. 
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"Redraw Map" refreshes the activity window. The map view may be un-zoomed back to the entire 
world by pressing "Reset Map", keeping the same mode and with no loss of data. 

11.3.11 Drawing Area Window Features. 

When the Map Operations dialog box pops up, a world map will be drawn in the Drawing Area 
Window. The map has several static features. The user should use these to aid his choice of domain 
parameters such as size and terrain data set. Lines are used for cartographic features. Coverage areas 
are shown as square panels and indícate what resolution of data is available for a given portion of the 
world. The dynamic features of the map are objects drawn by the user. These are features used to 
represent data input and other selections. This chart is a key to the map’s features. An example with 
most of the features is presented in Figure 48. 

Static Features 


Shape 

Style 

Weight 

Color 

Feature 

Line 

Solid 

Thick 

Black 

Equator 

Line 

Solid 

Thick 

Black 

Prime Meridian 

Line 

Solid 

Thin 

Black 

Tropics 

Line 

Dash 

Thin 

Black/Gray 

lOdeg Latitude 

Line 

Dash 

Thin 

Black/Gray 

lOdeg Longitude 

Space 

— 

— 

White 

Areas of 5 minute Resolution Coverage 

Square 

Filled 

Thin 

Yellow 

Areas of 30 second Resolution Coverage 



Figure 48. Draw area features. 


125 







Dynamic Features 
Shape Style 
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Feature 


Región using 5 minute Resolution Data 
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Refinement Disk 
ADM Source Location 


11.4 Set Time Domain. 

When performing a simulation of the atmosphere the user creates a weather forecast. Even if the user 
simulates an historie event, the user still starts with the atmosphere in some State and evolves it in time 
to an end State. The Set Time Domain dialog box is the mechanism used to set the temporal bounds on 
the simulation. 


The first operation in the Operations Column is "Set Time Domain". This placement reflects the idea 
that the time bounds of an OMEGA simulation are the most fundamental part. The user first decides the 
time period of the simulation before setting other parameters. When the button is pressed the Set Time 
Domain dialog box (Figure 49) pops up. 

This dialog box is one of the easiest dialogs to use. The dialog consists of three sections: Forecast Time 
Settings and two rows of Action Buttons. How to use each section is described, followed by exact 
specifications of each control’s actions, valid data range, and default valué. 

The data for each región is stored intemally as both the times you input and the total length of time for 
the simulation. 


11.4.1 Forecast Time Settings. 


This section has eight text entry boxes. By clicking the first mouse button within the box the user enters 
his choices. The cursor location is a blinking bar. The user may now type as usual in the box to enter 
the desired Start and End Times in broken-down form for the simulation. 

On most systems both backspace and delete will work to remove digits. Be sure to abide by the digit 
conventions Usted in parenthesis beside the box when entering data. 


Start Year - Year of Start Time in four digit form. 
Valid Range: 1900 to 2100 
Default Valué: current year 

Start Month - Month of Start Time in two digit form. 
Valid Range: 1 to 12 
Default Valué: current month 


Start Day - Day of Start Time in two digit form. 
Valid Range: 1 to 31 
Default Valué: current day 


126 




Start Hour - 


End Year - 


End Month - 


End Day - 


End Hour - 



Figure 49. Set time domain dialog box. 

Hour of Start Time in four digit form. 

Valid Range: 00 to 23 
Default Valué: 00 

Year of Start Time in four digit form. 

ValidRange: 1900 to 2100 
Default Valué: current year + 24 hours 

Month of Start Time in two digit form. 

Valid Range: 1 to 12 

Default Valué: current month + 24 hours 

Day of Start Time in two digit form. 

Valid Range: 1 to 31 

Default Valué: current day + 24 hours 

Hour of Start Time in four digit form. 

Valid Range: 00 to 23 
Default Valué: 00 
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11.4.2 Action Buttons. 


Clear Form - This button resets the text en tries to blanks. 

Default - This button resets all text entries to their default valúes. 

OK - This button checks the Time data, cause XOMEGA to write out a "time.inp" file in the case 

directory, and cióse the dialog box. If any valúes are in error, a Set Time Domain Error dialog 
box will pop up alerting the user of the specific errors and suggestions on how to correct 
them. See Section 11.4.3 for information on the Error dialog box. No "time.inp" file will be 
written if there are errors in the data. 

Cancel - This button declines the Time data, reset the valúes to default, and cióse the dialog box. 

There is no confirmation process. 

Help - This button opens the Help dialog and load this section. 

11.4.3 Set Time Domain Error. 

If an error is made while entering valúes into the text entries, a dialog box pops up alerting the user of 
the errors and suggesting ways to correct them. The top part of the dialog box contains the error 
message and the suggestions. The bottom part has a row with two action buttons: 

OK - This button dismisses the Error dialog. 

Help - This button will open the Help dialog and load this section. 

11.5 Set Spatial Domain. 

The domain selection process is one of the most critical steps in creating an OMEGA simulation. If the 
domain parameters are not set correctly the model quickly produces errors and may stop altogether. 
When "Set Spatial Domain" is pressed two dialog boxes. Set Spatial Domain and Map Operations, pops 
up and a world map is displayed in the Drawing Area window. The dialog boxes are used to input 
spatial domain information, which then displays on the map. 

The Set Spatial Domain dialog box (Figure 50) is the main control for the OMEGA domain. It is used 
to define the horizontal regions of the domain and to define the vertical spacing of the grid levels. The 
dialog has five sections: Current Región Controls, Región Parameters, Vertical Parameters, and two 
rows of Action Buttons. Each section is described below, followed by exact specifications of each 
control’s actions, valid data range, and default valué. 

The data for each región is stored intemally as a list of regions. The Current Región’s information is 
displayed in the Set Spatial Domain dialog. To view data from the other regions the Current Región 
Controls are used. 

Please see Chapter 5 for more information on the OMEGA Grid Generator. 

11.5.1 Current Región Controls. 

When using the cursor to endose a región on the map, the numbers located at the top of the Set Spatial 
Domain dialog change accordingly. By clicking on the arrow buttons the user may cycle through the 
regions. On the map the current región is outlined by a thick rectangle. Any other regions the user may 


128 



have defined will be drawn with thin rectangles. Any región may be deleted by clicking on "Delete" 

(Note: there is NO confirmation process). 

Current Región - A dynamic label showing the current región and the total number of regions. The 
current región’s valúes appear in the Región Parameters’ text entry boxes and radio 
buttons. 

Left Arrow - Decrement the Current Región filling the Región Parameters’ information and 
drawing the Current Región with a thick rectangle Cyclic action going down: 1 <- 2 <- 
3 <- 1 <- 2 

Right Arrow - Increment the Current Región filling the Región Parameters’ i n f ormation and drawing 
the Current Región with a thick rectangle. Cyclic action going up: 1 -> 2 -> 3 -> 1 -> 
2 ... 



Figure 50. Set spatial domain dialog 
box. 
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Delete - Remove the Current Región and shift up any subsequent Regions. For example: four regions 
are defined, the current región number 3, and Delete is clicked. Región 3 is removed from 
the list, and Región 4 becomes Región 3. 

11.5.2 Región Parameters. 


The Región Parameters section displays the bounding coordinates of the región in latitude and longitude, 
the desired resolution range in kilometers, and the data set used for that región. Each región will be 
displayed on the map by a colored rectangle. The color is determined by the terrain data set chosen. 

The user must ensure that the requested data set is available for the part of the world being simulated 
(e.g., a simulation of the mid-Pacific may only use the 5 minute (white area on the map) data set; 
simulations of the United States may use 30 second (yellow panels on the map) data). 

When making an OMEGA domain it is best to choose the outer-most región with 5-minute data and then 
make smaller sub-regions with the 30-second data. These regions must be nested together, and must 
have resolution ranges that slightly overlap, so that the final grid has a smooth transition from coarse to 
fine resolution. Note that no subsequent región may be drawn outside Región 1, so the first región must 
endose the entire intended domain. 

As an example, consider a simulation for Washington DC. A total domain (Región 1) covering the Mid- 
Atlantic part of the US might be specifíed, using the 5-minute data set with a resolution range of 50- 
70km. Then, Región 2 covering Maryland and Virginia using the 30-second data set with a resolution 
range of 35-55km and Región 3 covering the Chesapeake Bay area using the 30-second data set with a 
resolution range of 10-36km could be defined. Finally, Región 4 consisting of the city itself with 30- 
second data set with a resolution range of 3-1 Ikm could be specifíed. The final grid will be smooth with 
a reasonable spatial resolution gradient. 


Minimum Longitude - 
Máximum Longitude - 

Minimum Latitude - 


Máximum Latitude - 


Min. Resolution - 


Max. Resolution - 


Western edge of Current Región. Negative valúes are West Longitude. 
Valid Range: -180.00 to 180.00 
Default Valué: 0.00 

Eastem edge of Current Región. Negative valúes are West Longitude. 
Valid Range: -180.00 to 180.00 
Default Valué: 0.00 

Southern edge of Current Región. Negative valúes are South Latitude. 
Valid Range: -90.00 to 90.00 
Default Valué: 0.00 

Northern edge of Current Región. Negative valúes are South Latitude. 
Valid Range: -90.00 to 90.00 
Default Valué: 0.00 

Minimum acceptable resolution for Current Región in km. 

Valid Range: 1 to 200 km 
Default Valué: 30 km 

Máximum acceptable resolution for Current Región in km. 

Valid Range: 1 to 200 km 
Default Valué: 60 km 
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Terrain Data Resolution - Radio button pair: only one button may be active. This button Controls 

the terrain data source to build the part of the grid defined by this región. 
Regions using 5-minute data will be drawn in black. Regions using 30- 
second data will be drawn in blue. Default Valué: 5 minute 

Shoreline Data Resolution - Radio button pair: only one button may be active. This button Controls 

the shoreline data source used to build the part of the grid defined by this 
región. 

Default Valué: 5-minute 

11.5.3 Vertical Parameters. 

The Vertical Parameters section allows the user to change how many total levels to use in the model, 
and the base thickness of Level 1. From these numbers the total height of the domain is calculated and 
displayed on the dialog in kilometers. For most applications the default valúes should be used. 

Total Number of Levels in Domain - This is the total number of levels of the computational domain. 

Valid Range: 5 to 100 
Default Valué: 30 

Level 1 Thickness - This is the thickness in meters of the first level. 

Valid Range: 5.0 to 250.0 meters 
Default Valué: 30 

11.5.4 Action Buttons. 

Clear Form - This button resets the text entries to zeros, de-select terrain and shoreline data sets, and 

reset the Vertical Parameters to default valúes. It does not modify the Región Data in any 
way. The user must modify the valúes before continuing. 

Default - This button resets all text entries and radio buttons to their default valúes. It does not 

modify the Región Data in any way. The user must modify the valúes before continuing. 

OK - This button checks the Región data, cause XOMEGA to write out a "grid.inp" file in the 

grid directory and cióse the dialog box. If any valúes are in error, a Set Spatial Domain 
Error dialog box will pop up alerting the user of the specific errors and suggestions on 
how to correct them. See Section 11.5.5 for information on the Error dialog box. No 
"grid.inp" file will be written if there are errors in the data. 

Cancel - This button declines the Región data, remove all regions, reset the valúes to default, and 
cióse the dialog box. There is no confirmation process. 

Help - This button opens the Help dialog and load this section. 

11.5.5 Set Spatial Domain Error. 

While entering valúes into the text entries if the user makes errors a dialog box pops up alerting him of 
the errors and suggesting ways to correct them. The top part of the dialog box contains the error 
message and the suggestions. The bottom part has a row with three action buttons: 

Delete Current Región - This button will remove the Current Región from the list similar to Delete on 

the Set Spatial Domain dialog. Be sure that the Current Región is the one 
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containing the error before clicking this button. Use the arrow buttons if 
needed to change to the Current Región. This button will dismiss the Error 
dialog. 

Modify Regions - This button will retum control to the Set Spatial Domain dialog box and 

dismiss the Error dialog box. You may now modify the regions as usual. 

Help - This button opens the Help dialog and load this section. 

11.5.6 Set Refinement Dialog. 

XOMEGA allows for an advanced feature cali Point Refinement. This is activated using the 
“Advanced” menú and “Set Refinement Disks”, or the “Map Operations” dialog and “Refinement 
Disk”. The cursor is then used to draw circles. These circles represent regions of the domain whose 
resolution will be fixed to the specified valué. The Grid Generator makes sure that the cells produced 
within the disk are of a fixed resolution. 


Selecting a refinement disk enables the user to have increased resolution about a specific point, without 
defining another región. This saves on I/O process time, allowing the user to pinpoint a specific 
resolution instead of a range. 

The effect of point refinement is Gaussian in nature: cells produced cióse to the center point of the disk 
are as cióse to the prescribed resolution as possible, while those farther out will vary in resolution to 
match the bounding region's resolution. 

The Valid Range for the resolution of a Refinement Disk is determined by the bounding Región. The 
resolution must be no smaller than 0.5 times the Minimum Resolution of the bounding región, and no 
greater than the Máximum Resolution of the bounding región. For example if the user has defined a 
región with 10-20km resolution and wants to put a refinement disk within this región, he must make sure 
the resolution of the disk is between 5-20km. This ensures a smooth transition from the refinement disk 
to the bounding región. 

The Set Refinement Disks dialog (Figure 51) has four sections: Current Disk Controls, Disk 
Parameters, and two rows of Action Buttons. 


Current Disk - A dynamic label showing the current disk and the total number of disks. The current 
disk's valúes appear in the Disk Parameters' text entry boxes. 

Left Arrow - Decrement the Current Disk filling the Disk Parameters' information and drawing the 
Current Disk in thick circle. 

Cyclic action going down: 1 <- 2 <- 3 <- 4 <- 1 ... 

Right Arrow - Increment the Current Disk filling the Disk Parameters' information and drawing the 
Current Disk in thick circle. 

Cyclic action going up: 1 -> 2 -> 3 -> 4 -> 1 ... 


Delete - 


Remove the Current Disk and shift up any subsequent Disks. 
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Longitude - 

Latitude - 

Radius - 

Resolution - 

Clear Form - 
Default - 

OK- 


Curoent Disk 1 of 1 


Pelete 


Disk Parameters: 


Longitude •]—78.99 

* 

Latitude ÍÍ¡7~81~ 


Radius (km) |]53.18 
Resolution (km) 


Í3B.00 



Defauhs 


Clear Form 



Done 

Cancel 

ÍMp 



Figure 51. Set refinement disks dialog box. 

Longitude valué of the Current Disk’s center point. 

ValidRange: -180.00 to 180.00 
Default Valué: 0.00 

Latitude valué of the Current Disk’s center point. 

Valid Range: -90.00 to 90.00 
Default Valué: 0.00 

The radius of the Current Disk in kilometers. 

Valid Range: 0.5 km to 100 km 
Default Valué: 0.00 

The desired resolution for the Current Disk. 

Valid Range: See discussion above 
Default: 10.00 

This button resets the text entries to zeros. It does not modify the Disk Data in 
any way. The user must modify the valúes before continuing. 

This button resets all text entries and radio buttons to their default valúes. It does 
not modify the Región Data in any way. The user must modify the valúes before 
continuing. 

This button checks the Disk data, cause XOMEGA to write out a "grid.ref' file in 
the grid directory and cióse the dialog box. If any valúes are in error, a Set 
Refinement Disk Error dialog box will pop up alerting the user of the specific 
errors and suggestions on how to correct them. No "refine.lst" file is written if 
there are errors in the data. 











This button declines the Disk data, remove all disks, reset the valúes to default, 
and cióse the dialog box. There is no confirmation process. 

This button opens the Help dialog and load this section. 

11.6 Set ADM Sources. 

OMEGA has a built-in Atmospheric Dispersión Model (ADM) and XOMEGA allows the user to select 
source locations for this model. All these locations must be within the model domain. When "Set ADM 
Sources" is clicked two dialog boxes, Set ADM Sources and Map Operations, pops up and a world map 
is displayed in the Drawing Area window. The dialog boxes are used to input ADM source Information 
that will then be displayed on the map. 

The Set ADM Sources dialog box (Figure 52) is the main control for the ADM program. It is used to 
define the source locations and strengths. The dialog has five sections: Atmospheric Dispersión Model 
Options, Current Source Controls, Source Parameters, and two rows of Action Buttons. Each section 
will be described below, followed by exact specifications of each control’s actions, valid data range, and 
default valué. 

The data for each source is stored intemally as a list of sources. The Current Source’s information is 
displayed in the Set ADM Sources dialog. Data from the other sources is viewed using the Current 
Source Controls dialog. 

Note that several of the Controls are only used with advanced ADM modes. The general user is 
encouraged to leave the default valúes as is, and specify only Location, Altitude, and Time valúes. 

Please see Chapter 8 for more information on the Atmospheric Dispersión Model. 

11.6.1 Atmospheric Dispersión Model Options. 

This group of option menus is responsible for setting the most general features of the ADM module. 
These options will apply to all sources selected. The options control which ADM program is run, and 
how it behaves. 

Mode - Currently this option menú has two choices: Off and Tracer. Off will cause the ADM 
program to be skipped when running the OMEGA Model. Tracer actives tracer type 
particles in the ADM model. 

Dcfáult' Trsccr 

Diffusion - This option menú has two choices: Puff and LPDM. Puff causes the ADM program to use 
puff computations when it runs. LPDM causes the ADM program to use Lagrangian 
Particle Diffusion when it runs. 

Default: Puff 

Eulerian - This option menú has two choices: Off and On. If On is selected, The ADM takes the ñrst 
three sources defined and duplicates them as Eulerian Sources. 

Default: Off 


Cancel - 
Help - 
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Figure 52. Set ADM sources dialog box. 


11.6.2 Current Source Controls. 

When using the cursor to click on an ADM Point on the map, the numbers located in this section 
increments. Clicking on the arrow buttons allows the user to cycle through the sources. On the map the 
current sources are shown as an open orange dot. Any other sources the user may have defined are 
drawn with solid orange dots. The user may delete any source by clicking on "Delete" (Note: there is 
NO confirmation process). 

Current Source - A dynamic label showing the current source and the total number of sources. The 
current source’s valúes appear in the Source Parameters’ text entry boxes. 

Left Arrow - Decrement the Current Source filling the Source Parameters’ information and drawing 
the Current Source as an open orange dot. 

Cyclic action going down: 1 <- 2 <- 3 <- 4 <- 1 ... 
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Right Arrow - Increment the Current Source filling the Source Parameters’ information and drawing 
the Current Sources an open orange dot. 

Cyclic action going up: 1 -> 2 -> 3 -> 4 -> 1 ... 

Delete - Remove the Current Source and shift up any subsequent sources. For example: four 

sources are defined, the current source is number 3, and Delete is clicked. Source 3 is 
removed from the list, and source 4 becomes source 3. 

11.6.3 Source Parameters. 


The Source Parameters section displays all the information about a source, including its location, time 
bounds, and strength. Typical valúes for some parameters are listed for 1000-micron diameter dust 
particle. Other types of particles will have different typical valúes. 


Longitude - 
Longitude - 

Aititude - 

Start Time - 

End Time - 

Injection Interval - 

Diameter - 

Mass - 

Density - 

Vapor Concentradon - 


Longitude of current source location. 

ValidRange: -180.00 to 180.00 

Default Valué: 0.00 

Latitude of current source location. 

Valid Range: -180.00 to 180.00 
Default Valué: 0.00 

Aititude of current source location in meters above ground level. 

ValidRange: 1.00 to lOO.OOm 
Default Valué: 20.00 

Number of seconds after start of simulation to begin release of particles from 
current source. 

Valid Range: 0.00 to máximum number of seconds in simulation 
Default Valué: 0.00 

Number of seconds start of simulation to halt release of particles from current 
source. 

Valid Range: 0.00 to máximum number of seconds in simulation 
Default Valué: máximum number of seconds in simulation, if defined. 
Number of seconds between release of particle from current source. 

Valid Range: 1.00 to 10800.00 seconds (30 hours) 

Default Valué: 1800.00s (1/2 hour) 

Diameter of particle in meters released from current source. 

Default Valué: 0.00 m 
Typical Valué: 0.001 m 

Mass of particle in kilograms released from current source. 

Default Valué: 0.00 kg 
Typical Valué: 1.05e-6kg 

Density of particle in kilograms per cubic meter for current source. Note this 
number must be consistent with the Mass and Diameter. 

Default Valué: 0.00 kg/m3 
Typical Valué: 2000.00 kg/m3 

Vapor Concentration of puff in kilograms per cubic meter for current source. 
Note this number must be consistent with the Mass and Diameter. 

Default Valué: 0.00 kg/m3 
Typical Value: 1.00 kg/m3 
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Sigma in x - 


Initial sigma x-value of the puff in meters for the current source. 

Default Valué: 0.00 m 
Typical Valué: 5.00 m 

Sigma in y - Initial sigma-y valué of the puff in meters for the current source. 

Default Valué: 0.00 m 
Typical Valué: 5.00 m 

Sigma in z - Initial sigma-z valué of the puff in meters for the current source. 

Default Valué: 0.00 m 
Typical Valué: 5.00 m 

Particle Multiplier - Number of actual particles to release from current source. If this valué is not 

one, it is used to modify the source term in the ADM so that múltiple 
particles are released at each interval from the same source. For example if 
this valué is three, then the source releases three particles every interval, not 
just one. 

Valid Range: 1 to 10 
Default Valué: 1 

11.6.4 Action Buttons. 

Clear Form - This button resets the text entries to zeros, and re-set the option menus to the default 

valúes. It does not modify the Source Data in any way. The user must modify the valúes 
before continuing. 

Default - This button resets all text entries and option menus to their default valúes. It does not 

modify the Source Data in any way. The user must modify the valúes before continuing. 

OK - This button checks the Source data, causes XOMEGA to write out a "adm.inp" file in the 

case directory and cióse the dialog box. If any valúes are in error, a Set ADM Sources 
Error dialog box will pop up alerting the user of the specific errors and suggestions on 
how to correct them. See Section 11.6.5 for information on the Error dialog box. No 
"adm.inp" file is written if there are errors in the data. 

Cancel - This button declines the Source data, remove all sources, reset the valúes to default, and 
cióse the dialog box. There is no confirmation process. 

Help - This button opens the Help dialog and load this section. 

11.6.5 Set ADM Sources Error. 

If an error is made while entering valúes into the text entries, a dialog box pops up alerting the user of 

the errors and suggesting ways to correct them. The top part of the dialog box contains the error 

message and the suggestions. The bottom part has a row with three action buttons: 

Delete Current Source - This button removes the Current Source from the list similar to Delete on the 

Set ADM Sources dialog. Be sure that the Current Source is the one 
containing the error before clicking this button. Use the arrow buttons if 
needed to change to the Current Source. This button dismisses the Error 
dialog. 

Modify Regions - This button retums control to the Set ADM Sources dialog box and dismiss 

the Error dialog box. You may now modify the sources as usual. 

Help - This button opens the Help dialog and load this section. 
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11.7 Run Grid Generator. 


When the user is satisfied with the arrangement of Spatial Domain Regions, the ADM Source locations 
(and possible Refinement Disks), the next step is to create the OMEGA grid. The program that 
translates the 2-D domain región information into a full 3-D computational mesh is called the OMEGA 
Grid Generator. 

Since this program is independent of XOMEGA it needs to be called with a controlled process. The 
mechanism for doing this is to use a UNIX pipe to send the text output of the Grid Generator back to 
XOMEGA into a read-only text window. While the Grid Generator is running all of its text output is 
captured and displayed so that the user may watch the progress. 

The Grid Generator first reads the file "grid.inp" to get information about the domain regions and levels. 
The grid is then made by first gathering the requested data and "stitching" together the individual data 
sets into a seamless block of terrain and shoreline data. Next, a series of triangles is interpolated to the 
underlying data. The interpolation process adds and deletes triangle vértices until all triangles are within 
the resolution range requested. The grid is now a single entity with varying resolution, but no distinct 
regions. The regions are temporary boundaries used to specify desired resolution ranges and data set 
usage. Please see Chapter 5 for more information on the OMEGA Grid Generator. 

When the user clicks on "Run Grid Generator" a large dialog box pops up. This dialog consists of two 
sections: Text Output Window and a row of Action Buttons. Figure 53 shows the dialog after the Grid 
Generator has started. 

Note that if the Spatial Domain has not been specified the Run Grid Generator dialog will not pop up. 
Instead an error message will pop up informing the user of the error. When acknowledging the error by 
clicking the "OK" button the Set Spatial Domain dialog will pop up. Please see Section 11.5 for 
information on setting the Spatial Domain. 

11.7.1 Text Output Window. 

When initially popped up the Text Output Window is blank. As the grid generation progresses the 
window becomes filled with text. The output is saved as it scrolls off the window, and even saved when 
the dialog is popped down. The scroll bar on the right hand side of the window may be used to view 
text beyond the window. On subsequent pop-ups the previous information is still in the Text Output 
Window. Any new text will be appended in the window. 

11.7.2 Action Buttons. 

Cióse - This button dismisses the dialog box. XOMEGA will read the report file created 

by the Grid Generator, "grid/grid.rep", and pop up a dialog with this information 

Compression Ratio: Ratio of simulated time to cpu time for the grid 24-Hours Run Time: Real time to 
simúlate 24 hours with this grid. Note that this information is machine dependent. 

If the Grid Generator is unable to make these estimations, the dialog box contains the valúes "-99.0" for 
the numbers. Please see Chapter 5 for more information on the OMEGA Grid Generator. 
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********* BEGINNING OMEGA GRID GENERATION ********* 


ARRAY ALLOCATIONS: 
MVM: 2508 MEM: 

MEMORY REQUIRED: 
SETUP: Requested 


5016 MRM: 31 


5630881 words. 


file-> grid,ref, exists 

grid-generation will inelude pt. refinement. 

bb + b^-k-k-k-k'h'irii'k-k'k-k-k-kJr-k'k’b'k-k’k'k-k-k-k’ir’k-k-k'k-k-kir-k-k-ir-k-ir-jc-iT-k-k-je-k-k-lc-k-ie-k-ir-k-k-k-k’k’k-k-k 

* Grid generation will be based on cell areas rather than * 

* edge lengths. * 

* Celis having small aspect raties may have edge lengths * 

* outside the specified resolution range. * 


¡ PROCESSING SMIN TERRAIN DATA, DOMAIN 1 
| /nf s/tornado/orne gal/orne ga/te r r ain/te r5min/te r5min. dat 

| PROCESSING SMIN LAND/WATER DATA, DOMAIN 1 

| 

1 DOMAIN 2 BAND 1 

\ reading 30sec terrain panel -80 35 

] OPENING 30SEC TERRAIN DATA FILE: n03Sw080. mvk 


Cióse 



Start - 


Figure 53. Run grid generator dialog box. 

When this button is pressed a check is made to determine if a grid file already 
exists. If a valid grid already exists the Grid Generation waming dialog box pops 
up. See Section 11.7.3 for information on the Waming dialog box. 


If no grid exists then the Grid Generator is started. While the Grid Generator is running only the Stop 
button will accept input. All other Controls on this dialog, and XOMEGA are locked out. 


Stop - 


If the Grid Generator is not running this button has no effect. If the Grid 
Generator is running this button will stop it. 


Note that when this button is pressed its effeets may not be immediately apparent. It may take a few 
seconds for the "Stop" command to propágate to the Grid Generator program, and for XOMEGA to reset 
its Windows and Controls. 


Since the extemal program is run through a UNIX pipe its output is sent back down the pipe to the 
dialog box. When the dialog is waiting for output from the program XOMEGA is unable to update its 
Windows or process a "Stop" button press. If Windows are covered and uncovered, parts of XOMEGA 
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may become gray. When the next line of output from the extemal program is piped to the dialog the 
Windows are redrawn and the button press processed. XOMEGA resets itself. 

Help - This button will open the Help dialog and load this section. 

11.7.3 Grid Generation Warning. 

This dialog box is popped up when the user tries to run the Grid Generator and a valid grid file already 
exists. The top part of the dialog contains the warning message. The bottom part has a row of three 
action buttons: 

Overwrite - This button starts the grid generation process. All existing information in the "grid/" sub- 
directory will be overwritten as the new grid is made. The output from the Grid Generator 
appears in the Text Output Window normally. If the program has no errors the Run Grid 
Generator operation is completed. The Grid Generation Warning dialog is dismissed. The 
Run Grid Generator dialog remains. 

Proceed - This button accepts the existing grid. All existing grid information is left untouched and 
the Run Grid Generator operation is completed. The Grid Generation Warning dialog is 
dismissed. The Run Grid Generator dialog is dismissed. 

Cancel - This button declines the existing grid. All existing grid information is left untouched and 
the Run Grid Generator operation is NOT completed. The Grid Generation Warning dialog 
is dismissed. The Run Grid Generator dialog is dismissed. 

11.8 Set Meteorology. 

An OMEGA simulation requires weather data to build initial conditions and boundary conditions. 
Common weather data sources are models whose domains are global or hemispheric, and whose forecast 
times are days. Observational data may also be used. Both surface and upper air observations may be 
used in an OMEGA run. 

When constructing an OMEGA simulation the user must be sure that the meteorology data is valid for 
the time and space of your simulation. Not all models have global coverage and some model’s forecast 
times may not match the user’s times. 

XOMEGA uses the concept of meteorology segments to allow for máximum flexibility when specifying 
the input meteorology. Essentially, a meteorology segment is a subdivisión of the time domain coupled 
with the sources of data used for that time period. A meteorology segment is defined by parameters that 
allow the specification of different sources of meteorology data be used for different times. 

As an example, consider modeling the weather for the city of Washington DC. The user might want to 
produce a 3-day forecast for the Metro región. Initial conditions for 00Z on the first day will have to be 
generated as well as boundary conditions for all days of the simulation. 

This city is within the domain of nearly all major atmospheric models, so grid model data should be 
available for the user’s domain. Suppose the user is interested in small-scale effeets and thus wants to 
use the highest resolution data possible to generate the boundary conditions. Most high-resolution 
models do not forecast for total forecast period (3-days) so the same model cannot be used for all of the 
boundary conditions. 
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This is where meteorology segments become very useful. The user can defíne a series of meteorology 
segments that specify which data sources to use for different periods of his forecast. Meteorology data 
can be mixed and matched to use the highest resolution possible at each analysis time. 

Assuming he has both a NWS feed and a FNMOC feed, his forecast period may be split into these three 
segments: 

Segment 1 covers the first 24-hour period and use ETAXP data. 

Segment 2 covers the next 24-hour period and use NORAPSN data. 

Segment 3 covers the last 24-hour period and use NOGAPS data. 

Using meteorology segments gives the user the flexibility to generate boundary conditions from 
different data sources at different times. This allows the use of the best data available in OMEGA 
simulations. The Set Meteorology dialog box (Figure 59) is the control for defining meteorology 
segments. The dialog has five sections: Time Domain label, Current Segment Controls, Segment 
Parameters, and two rows of Action Buttons. How to use each section will be described, followed by 
exact specifications of each control’s actions, valid data range, and default valué. 

Note that if the Time Domain data has not been filled the Set Meteorology dialog does not pop up. 
Instead an error message will pop up informing you of the error. When acknowledging the error by 
clicking the "OK" button the Set Time Domain dialog will pop up. Please see Section 11.4 for 
information on setting the Time Domain. 

11.8.1 Time Domain Label. 


This section of the dialog is a simple reminder of the Time Domain. It is presented so that the user has 
the total time span of your simulation available as you define meteorology segments. If the Time 
Domain is changed, the labels are updated. 

11.8.2 Current Segment Controls. 


Initially the dialog has only one segment defined spanning the entire time domain. The Controls in this 
section are used to add, delete and cycle through all defined regions. 


Current Segment - A dynamic label showing the current segment and the total number of segments. 

The current segment has its Segment Number Label non-shaded. All other Segment 
Numbers are shaded. 


Left Arrow - 
Right Arrow - 
Add - 


Delete - 


Decrement the Current Segment. 

Cyclic action going down: 1 <- 2 <- 3 <- 4 <- 1 ... 

Increment the Current Segment. 

Cyclic action going up: 1 -> 2 -> 3 -> 4 -> 1 ... 

Add a segment to the end of the list. The text entry and option menus will be set to 
default valúes. It is the user’s responsibility to adjust the Time entries to fit within 
the bounds. Up to five meteorology segments may be defined. 

Remove the Current Segment and shift up any subsequent Segments. If four 
segments are defined, and Segment 3 is current, Delete removes Segment 3 from 
the list, and Segment 4 becomes Segment 3. 
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Figure 54. Set meteorology dialog box. 


11.8.3 Segment Parameters. 


Segment Number - 


Start Time - 


End Time - 


Gridded Weather Model - 

Surface Observations - 

Rawinsonde Observations - 


Dynamic label used to highlight which segment is the current segment. 

The current segment will have its number in bold print. All other 
segments will have shaded gray numbers. 

Valid Range: 1 to 5, in fíxed sequence 
Default Valué: 1 

Start Time of Segment in full digit form: four digits for year, two digits 
for month, day and hour, single space between numbers. 

Valid Range: 1900 01 01 00 to 2100 12 31 23 

Default Valué: valué copied from Start Time of Time Domain 

End Time of Segment in full digit form: four digits for year, two digits for 

month, day and hour, single space between numbers. 

Valid Range: 1900 01 01 00 to 2100 12 31 23 

Default Valué: valué copied from End Time of Time Domain 

Source model to use for Segment. 

Valid Range: Installation dependent 
Default: Installation dependent 
Source of observations to use for Segment. 

Valid Range: Installation dependent 
Default: None 

Source of upper air observations for Segment. 

Valid Range: Installation dependent 
Default: None 
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11.8.4 Action Buttons. 


Clear Form - This button deletes all segments, reset the text entries to blank and de-select all option 
menus. 

Default - This button deletes all segments, except the first. It will set the first segment to default 
valúes. 

OK - This button checks the Segment data, cause XOMEGA to write out a "met.inp" file in 

the prep directory and cióse the dialog box. If any valúes are in error, a Set Meteorology 
Error dialog box will pop up alerting the user of the specific errors and suggestions on 
how to correct them. See Section 11.8.4 for Information on the Error dialog box. No 
"met.inp" file will be written if there are errors in the data. 

Cancel - This button declines the Segment data, remove all segments, clear all valúes, and cióse 

the dialog box. There is no confirmation process. 

Help - This button opens the Help dialog and load this section. 

11.8.5 Set Meteorology Error. 

If an error is made while entering valúes into the text entries, a dialog box pops up alerting the user of 
the errors and suggesting ways to correct them. The top part of the dialog box contains the error 
message and the suggestions. The bottom part has a row with two action buttons: 

OK - This button acknowledges the error and dismiss the Error dialog. 

Help - This button opens the Help dialog and load this section. 

11.9 Fetch Meteorology. 

When the user is satisfied with the subdivisión of his time domain into one or more meteorology 
segments, the next step is to gather together the raw data. The program that searches the file system for 
the appropriate meteorology data files is called Fetch. 

Since this program is independent of XOMEGA it needs to be called with a controlled process. The 
mechanism for doing this is to use a UNIX pipe to send the text output of Fetch back to XOMEGA into 
a read-only text window. While the Fetch program is running all of its text output is captured and 
displayed so that the user may watch the progress. 

To keep the individual meteorology segments sepárate each segment is processed in its own 
subdirectory. These directories are below the prep directory. For example "case/prep/segl" will hold all 
the raw data specified in Segment 1. XOMEGA writes a simple shell script that creates the sub¬ 
directory and calis Fetch from there for each segment. Fetch is called by the script with command line 
arguments specifying the time span to use and what source of data to gather. Fetch scans the file system 
for the files. If the files are found a temporary copy is placed in the segment directory. If a file is 
missing, Fetch will look for another file that satisfies the time of the missing file. 

When the user cíicks on "Fetch Meteorology" a large dialog box pops up. This dialog consists of two 
sections: Text Output Window and a row of Action Buttons. Figure 55 shows the dialog after Fetch has 
started. 
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Figure 55. Fetch meteorology dialog box. 

Note that if the Meteorology Segments data has not been filled the Fetch Meteorology dialog will not 
pop up. Instead an error message will pop up informing the user of the error. When acknowledging the 
error by clicking the "OK" button the Set Meteorology dialog will pop up. Please see Section 11.8 for 
information on setting the Meteorology. 

11.9.1 Text Output Window. 

When initially popped up the Text Output Window will be blank. As the gathering of meteorology data 
progresses the window will become filled with text. The output is saved as it scrolls off the window, 
and even saved when the dialog is popped down. The scroll bar on the right hand side of the window can 
be used to view text beyond the window. On subsequent pop ups the previous information will still be 
in the Text Output Window. Any new text will be appended in the window. 

11.9.2 Action Buttons. 

Cióse - This button dismisses the dialog box. 

Start - When this button is pressed a check is made to determine if a set of meteorology data already 

exists for the case. If a valid set already exists the Fetch Meteorology Waming dialog box pops 
up. 

See Section 11.9.3 for information on the Waming dialog box. 

If no meteorology exists for the case then Fetch is started. While Fetch is running only the Stop button 
will accept input. All other Controls on this dialog, and XOMEGA are locked out. 

Stop - If Fetch is not running this button has no effect. 

If Fetch is running this button stops it. 
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Note that when you press this button its effects may not be immediately apparemt. It may take a few 
seconds for the "Stop" command to propágate to Fetch program, and for XOMEGA to reset its Windows 
and Controls. 


Since the extemal program is run through a UNIX pipe its output is sent back down the pipe to the 
dialog box. When the dialog is waiting for output from the program XOMEGA is unable to update its 
Windows or process a "Stop" button press. If the Windows are moved (or covered and uncovered) parts 
of XOMEGA become gray. When the next line of output from the extemal program is piped to the 
dialog the Windows will be redrawn and the button press processed. XOMEGA will reset itself. 

Help - This button opens the Help dialog and load this section. 

11.9.3 Fetch Meteorology Warning. 

This dialog box is popped up when the user tries to run Fetch and a valid set of meteorology data already 
exists for the case. The top part of the dialog contains the warning message. The bottom part has a row 
of three action buttons: 


Overwrite - 


Proceed - 


Cancel - 


This button starts the Fetch process. All existing information in the "prep/seg*" sub- 
directories will be overwritten as the new set of data is gathered. The output from Fetch 
appears in the Text Output Window normally. If the program has no errors the Fetch 
Meteorology operation is completed. The Fetch Meteorology Warning dialog is dismissed. 
The Fetch Meteorology dialog will remain. 

This button accepts the existing set of meteorology data. All existing meteorology data is 
left untouched and the Fetch Meteorology operation is completed. The Fetch Meteorology 
Warning dialog is dismissed. The Fetch Meteorology dialog remains. 

This button declines the existing set of meteorology data. All existing meteorology data is 
left untouched and the Fetch Meteorology operation is NOT completed. The Fetch 
Meteorology Warning dialog is dismissed. The Fetch Meteorology dialog remains. 


11.10 Run Preprocessor. 


After all raw meteorology data has been gathered for the meteorology segments, the next step is to 
process this data into initial conditions and boundary conditions. There are several programs that work 
on the data and collectiveíy they are known as the OMEGA Preprocessor. A parent script program has 
been written that binds the other programs together. Since all processing is performed by the programs 
the parent script calis it will be referred to as the Preprocessor. 


Please see Chapters 6 and 7 for more information on the OMEGA Preprocessor. 

Since this program is independent of XOMEGA it needs to be called with a controlled process. The 
mechanism for doing this is to use a UNIX pipe to send the text output of the Preprocessor back to 
XOMEGA into a read-only text window. While the Preprocessor is running all of its text output is 
captured and displayed so that the user may watch the progress. 


XOMEGA writes a simple script ("prep/prep.csh") to manage the processing of the meteorology 
segments and to keep the boundary conditions produced in time sequence. This script calis the 
Preprocessor once for each meteorology segment. Within each segment boundary conditions are made. 
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After all the segments are processed the scripí will rename the boundary files to form a sequence and 
copy them up to the case directory. 

When the user clicks on "Run Preprocessor" a large dialog box pops up. This dialog consists of two 
sections: Text Output Window and a row of Action Buttons. Figure 56 shows the dialog after the 
Preprocessor has started. 

Note that if the Fetch Meteorology operation has not been completed (i.e. no meteorology data was 
gathered) the Run Preprocessor dialog will not pop up. Instead an error message will pop up informing 
the user of the error. When acknowledging the error by clicking the "OK" button the Fetch Meteorology 
dialog will pop up. Please see Section 11.9 for information on Fetch Meteorology. 

11.10.1 Text Output Window. 

When initially popped up the Text Output Window will be blank. As the preprocessing progresses the 
window will become filled with text. The output is saved as it scrolls off the window, and even saved 
when the dialog is popped down. The scroll bar on the right hand side of the window may be used to 
view text beyond the window. On subsequent pop ups the previous information will still be in the Text 
Output Window. Any new text will be appended in the window. 

11.10.2 Action Buttons. 

Cióse - This button dismisses the dialog box. 

Start - When this button is pressed a check is made to determine if valid set of initial and boundary 

conditions already exist. If a valid set already exists the Run Preprocessor Waming dialog box 
pops up. See Section 11.10.3 for information on the Warning dialog box. 



j Start of prep.csh 


j Start of segment 1 

Ü 

• ]2;monsoon — /scratch/turner/aaaa/prep/segl PREPDAT will. execute 
' GENINIT will ejecute. 

: BOUNDARY CONDITIONS will be prepared 
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? I am now processing the OBSERVATIONAL data 
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Figure 56. Run preprocessor dialog box. 
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If no set exists then the Preprocessor is started. While the Preprocessor is running only the Stop button 
will accept input. All other Controls on this dialog, and XOMEGA are locked out. 

Stop - If the Preprocessor is not running this button has no effect. If the Preprocessor is running this 
button stops it. 

Note that when this button is pressed its effects may not be immediately apparent. It may take a few 
seconds for the "Stop" command to propágate to the Preprocessor program, and for XOMEGA to reset 
its Windows and Controls. 

Since the extemal program is run through a UNIX pipe its output is sent back down the pipe to the 
dialog box. When the dialog is waiting for output ffom the program XOMEGA is unable to update its 
Windows or process a "Stop" button press. If the user moves Windows (or covers and uncovers them) 
parts of XOMEGA become gray. When the next line of output from the extemal program is piped to the 
dialog the Windows will be redrawn and the button press processed. XOMEGA will reset itself. 

Help - This button opens the Help dialog and load this section. 

11.10.3 Run Preprocessor Warning. 

This dialog box is popped up when the user tries to run the Preprocessor and a valid set of initial and 
boundary conditions already exists. The top part of the dialog contains the warning message. The 
bottom part has a row of three action buttons: 

Overwrite - This button starts the Preprocessor. All existing information in the "prep/seg*" sub- 
directories will be overwritten as the new files are made. The output from the 
Preprocessor appears in the Text Output Window normally. If the program has no errors 
the Run Preprocessor operation is completed. 

The Run Preprocessor Warning dialog is dismissed. 

The Run Preprocessor dialog remains. 

Proceed - This button accepts the existing files. The Run Preprocessor operation is completed. 

The Run Preprocessor Warning dialog is dismissed. 

Cancel - This button declines the existing files. The Run Preprocessor operation is NOT 
completed. 

The Run Preprocessor Warning dialog is dismissed. 

The Run Preprocessor dialog is dismissed 

11.11 Run OMEGA Model. 

Once all preliminary setup steps are complete the user is ready to run the OMEGA Model. The Model is 
a single program that will use the grid the user has created and the preprocessed meteorology data 
gathered to simúlate the atmosphere. 

Please see Chapters 4,7,8 for more information on the OMEGA Model. 

Since this program is independent of XOMEGA it needs to be called with a controlled process. Also 
since the Model may take a long time to mn the process should be placed in the background. The 
mechanism for doing both these things is to use a small script program that will run the Model in the 
background and send status messages back to XOMEGA. As with the other extemal programs, a UNIX 
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pipe is used to send the text output of the script (status messages) back to XOMEGA into a read-only 
text window. The Model output (iteration information) is sent to a file called "model.log" in the case 
directory. 

When clicking on "Run Model" a small dialog box pops up. This dialog consists of two sections: Text 
Output Window and a row of Action Buttons. Figure 57 shows the dialog after the Model has started. 

Note that if the Run Grid Generator or Run Preprocessor operations have not been completed (i.,e. no 
grid was made or data processed) then the Run Model dialog do not pop up. Instead an error message 
will pop up informing the user of the error. When acknowledging the error by clicking the “OK” button 
the appropriate dialogs pop up. 

11.11.1 Text Output Window. 

When initially popped up the Text Output Window will be blank. As the Model starts up the window 
will become filled with text. The output is saved as it scrolls off the window, and even saved when the 
dialog is popped down. The scroll bar on the right hand side of the window may be used to view text 
beyond the window. On subsequent pop-ups the previous information will still be in the Text Output 
Window. When the Model is fully running the window will show the current status of the run. About 
every 5 seconds the information will be updated. 



Figure 57. Run model dialog box. 


148 


















11.11.2 Action Buttons. 


Cióse - This button dismisses the dialog box. 

Start - When this button is pressed a check is made to determine if valid output already exist. If a 
valid set of output files already exists the Run Model Waming dialog box pops up. 

See Section 11.11.3 for information on the Waming dialog box. 

If no set exists then the Model is started. While the Model is running only the Stop button will accept 
input. All other Controls on this dialog, and XOMEGA are locked out. 

Stop - If the Model is not running this button has no effect. If the Model is running this button will stop 
the status reporting process. The Model itself will continué to run in the background and must be 
stopped from a terminal with the kill command. 

Note that when you press this button its effects may not be immediately apparent. It may take a few 
seconds for the "Stop" command to propágate to the Model program, and for XOMEGA to reset its 
Windows and Controls. 

Since the extemal program is mn through a UNIX pipe its output is sent back down the pipe to the 
dialog box. When the dialog is waiting for output from the program XOMEGA is unable to update its 
Windows or process a "Stop" button press. If the Windows are moved (or covered and uncovered) parts 
of XOMEGA become gray. When the next line of output from the extemal program is piped to the 
dialog the Windows will be redrawn and the button press processed. XOMEGA will reset itself. 

Help - This button opens the Help dialog and load this section. 

11.11.3 Run OMEGA Model Warning. 

This dialog box is popped up when you try to mn the Model and a valid set of output files already exists. 
The top part of the dialog contains the waming message. The bottom part has a row of three action 
buttons: 

Overwrite - This button starts the Model. All existing output files will be overwritten as the new files 
are made. The output from the Model will appear in the Text Output Window normally. 
The Run Model Waming dialog is dismissed. The Run Model dialog remains. 

Proceed - This button accepts the existing files. The Run Model operation is completed. The Run 
Model Waming dialog is dismissed. The Run Model dialog is dismissed. 

Cancel - This button declines the existing files. The Run Model operation is NOT completed. The 
Run Model Waming dialog is dismissed. The Run Model dialog is dismissed. 
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Section 12 


OMEGA Output and Visualization 


12.1 OMEGA Output. 

The OMEGA Model produces a variety of output files depending on the mode of operation. The 
frequency of output is controlled by the user and generally occurs for each hour of simulation time. The 
standard meteorological files are suffixed with ".pkb", which stands for "packed binary". Similar files 
which contain the same data but at full Computer accuracy are generated at user defined intervals and 
used for restarting the model. These files are suffixed with ".out". If the user is running the ADM, a 
suite of other files may be generated as well. These files are currently written in ASCII. 

Binary format files are smaller and more efficient to read. Thus, it is our aim to convert the ASCII 
format files to a binary format in the future. However, there is a cross-platform issue that has meant 
frequent post-processing of these binary files when porting to different platforms. We have resolved this 
problem by implementing a new type of format called "packed binary". 

The packed binary format is a form of compression which results in files which are written and read on 
the byte level by routines written in the C programming language. These files are now completely 
portable and may be cross ported to any machine, including PCs. There is a small penalty in runtime, 
which is far made up for by savings in post-processing. 

There is an additional cost in accuracy which is irrelevant for visualization. In the past, data has been 
written using eight bytes to represent one number. With packed binary, this same number may be 
represented by as little as one byte. Accuracy diminishes as fewer bytes are used. In practice we use 
two bytes, which results in accuracy to five significant digits, which is sufficient for most purposes. 

This will require one quarter the amount of disk space used in the past. 

There is a requirement for files with full accuracy in order to restart the model. These files will continué 
to be produced as in the past in platform dependent binary, but less frequently than the routine files. 

In the future we will be converting the remainder of OMEGA generated files to packed binary format 
which will result in a much smaller database and less time required for post-processing. 

12.2 XGRID - OMEGA Visualization Tool. 

Due to the unstructured nature of the OMEGA grid, it is not easy to generate graphics with many of the 
common visualization packages. Henee SAIC designed and built a visualization tool for OMEGA to 
handle it’s unique characteristics. This has also allowed us to make continuous improvements and add 
new features to the display system with relative ease. 

XGRID is comprised of a graphical user interface (GUI), a viewing area for displaying the requested 
graphics, and a diagnostic text window. The first step in viewing OMEGA data is to open a .pkb file. 
This is done using a file selection browser launched from the File menú. XGRID will then search for 
and load its corresponding grid file and load the omega data into memory. At this point XGRID 
displays the grid terrain in the first of four possible modes: Layer - a horizontal section of the domain. 
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Additionally, the Field menú now becomes visible, allowing the user to select any of the meteorological 
fields, which have been generated from the model. A number of additional fields may be available, 
including some which are automatically calculated by XGRID upon loading the data; wind speed and 
direction, and mean sea level pressure. Surface characteristics will be available if the corresponding .srf 
file exists (Figure 58). 

In layer mode the user may select a variety of different options: using the mouse and clicking the cursor 
on different cells in the viewing area will display all information for that cell in the diagnostic window. 
Clicking on the vertical scale bar to the right of the viewing area will change the horizontal level which 
is displayed. The user may use the View menú to have the data contoured, overlay map information, or 
overlay wind barbs or arrows. The viewing area may be zoomed in or out, and text, boxes, title or dock 
can be displayed, some of which can be user controlled through sepárate files (Figures 59, 60). 



Figure 58. XGRID terrain view upon initial data load. 
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The user may view data as a vertical cross section. When this mode is selected, the user may click on 
the display and drag a line across the area of interest. When the button is released, XGRID will display 
a vertical cross section of the current field corresponding to this line (Figure 61). 



Figure 61. XGRID vertical cross section. 

Vertical point profiles may be generated for any field with vertical extent (Figure 62). 

Lastly, one may view a thermodynamic diagram called a SkewT-logP, which depicts a vertical profile of 
temperature and dewpoint (Figure 63). 

OMEGA ADM files may be displayed on XGRID showing the location of each particle, in layer or 
cross section mode (see section 10). Particle trajectories and field of regard may also be displayed 
(Figure 64). When the model is run in Puff mode, concentrations of Lagrangian tracers may be viewed 
(Figure 65). 

The new versión of XGRID has undergone extensive intemal change. Most memory allocation is now 
done dynamically, which negates the requirement for recompilation when viewing larger grids. This 
and the removal of all FORTRAN code have reduced the program memory requirements from 135 MB 
to about 2 MB. Extensive streamlining of the code will make it much easier to adapt in the future. 
Planned capabilities inelude reading structured grids from other models, animation on the fly, a new 
GUI with múltiple display Windows, and grid editing. Additionally, the new packed binary files will 
allow us to have only currently required data in memory. This will allow even greater speed and 
memory savings. 
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Figure 62. XGRID profile. 

12.3 OMEGA Simulation Image Maker (OSIM). 

While XGRID has the ability to create graphics guided by the user, there is a requirement for System 
generated graphics. XGRID has the capability of reading and executing a script file, but must do so 
under an actively running XGRID process with the GUI displayed. A new command line versión of 
XGRID has been created which reads a script file for its instructions. SAIC has created a web page on a 
remóte system using this program which will substantially reduce the required bandwidth for viewing 
desired data. 

The first step in the process is to access the OSIM web page (Figure 66). This page displays the runs 
and valid times available in the database, which the user must select from. 
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Figure 63. XGRID skewT-logP diagram. 
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Figure 64. Field of regard and particle trajectories. 



Figure 65. Concentration valúes computed from the ADM Lagrangian tracers. 
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Figure 66. OSIM stage one. 
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Figure 67. OSIM stage two. 







The second step is the selection of desired fields and display options (Figure 67). When this step is 
completed, the web server calis a program which generates an XGRID script, then calis the command 
line versión of XGRID. The requested graphics are then linked back to the web server, which displays 
them on the page (Figure 68). 

This new capability will drastically reduce the bandwidth required to get model results to users in the 
field. It is no longer necessary to transport large data files across the Internet in order to view a 
relatively small portion of that data. 


OMEGA - indi - 1998051100 - 199805111200 



Figure 68. OSIM stage three. 
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Section 13 


Conclusions 


OMEGA has grown tremendously over the last three years. Major changes have been made throughout 
the OMEGA modeling System, from the GUI to the Post-processors. The OMEGA GUI, XOMEGA, has 
been made more powerful, while maintaining its user-friendliness. We have added case management, 
greatly expanded the options and capabilities for specifying the meteorological data to be used for initial 
and boundary conditions, increased the possibilities for grid refinement through the interface, and added 
additional Controls on the ADM model. The grid generator was made more capable, adding point 
refinement routines and automating the generation of surface description fields (e.g., land use, soil type, 
vegetation). The OMEGA data pre-processor was improved adding more sources of data that can be 
accepted and adding options for mixing and matching data sources to create the initial and lateral 
boundary conditions. 

The OMEGA model itself was greatly improved. A new time-split solver was implemented so that the 
entire model would not be restricted by the smallest cell. This has greatly reduced run-time and 
improved performance. Dynamic grid adaptation was implemented so that the resolution can be 
maintained in a limited región of intense interest. And a first prototype parallel versión using MPI has 
been tested. In addition, major changes were made in the model physics including the PBL model and 
the inclusión of inhomogeneous surface properties, explicit cloud microphysics, statistical output 
routines to support the HPAC model with direct valúes of the variance of the OMEGA wind parameters 
over an extended (e.g., 1 hour) sample period, and updating the re-start capability to inelude all of these 
changes. The Atmospheric Dispersión Model (ADM) was also upgraded to inelude a probabilistic 
transport model based on the concentration variance. 

The model output format was changed to a machine independent, packed binary data format that is both 
more efficient in its storage requirement and considerably speeds up the post-processing of OMEGA 
output. The post-processing now ineludes quantitative verification routines (see Volume II) and a new, 
Web-based graphics visualization system was prototyped, the OMEGA Web Links (OWL) and the 
OMEGA Simulation Image Maker (OSIM). These systems permit custom product generation and 
display over the World Wide Web with very modest hardware; the prototype system used a 233 MHz 
Pentium II PC. 

OMEGA has again re-defined the state-of-the art in numerical weather prediction and atmospheric 
dispersión. It has been subjected to a suite of operational and field tests, which is the subject of Volume 
13 of this report. 
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4300 SAN MATEO BLVD, NE, SUITE A220 
ALBUQUERQUE, NM 87110-1260 
ATTN: R. NEWELL 

ITT INDUSTRIES 
ITT SYSTEMS CORPORATION 
ATTN: AODTRA/DASIAC 
1680 TEXAS STREET, SE 
KIRTLAND AFB, NM 87117-5669 
ATTN: DASIAC 
ATTN: DASIAC/DARE 


MESO INC. 

185 JORDAN ROAD 
TRO Y, NY 12180 
ATTN: JOHN ZACK 
ATTN: K.WAIGHT 
ATTN: P. E. PRICE 
ATTN: S.YOUNG 


MISSION RESEARCH CORPORATION 
ASTER DIVISION 

2629 REDWING DRIVE, SUITE 310 
FORT COLLINS, CO 80526 
ATTN: C. TREMBACK 
ATTN: M. WEISSBLUTH 


MISSION RESEARCH CORPORATION 
P. O. BOX 1257 
HUNTSVILLE, AL 35807 
ATTN: B. BAUER 

PACIFIC-SIERRA RESEARCH CORPORATION 
2901 28TH STREET 
SANTA MONICA, CA 90405-2938 
ATTN: H. BRODE 


SCIENCE APPLICATIONS INT'L CORPORATION 
2111 EISENHOWER AVENUE, SUITE 303 
ALEXANDRIA, VA 22314 
ATTN: J. COCKAYNE 

SCIENCE APPLICATIONS INTL CORP 
10260 CAMPUS POINT DRIVE 
SAN DIEGO, CA 92121-1578 
ATTN: T. JARRETT 


SCIENCE APPLICATIONS INTL CORPORATION 
1100 FIRST AVENUE, SUITE 300 
KING OF PRUSSIA, PA 19406 
ATTN: J. SONTOWSKI 

SCIENCE APPLICATIONS INT'L CORPORATION 
2109 AIR PARK ROAD, SE 
ALBUQUERQUE, NM 87106 
ATTN: J. MANSHIP 


LOGICON - RDA 

6940 S. KINGS HIGHWAY, SUITE 202 
ALEXANDRIA, VA 22310 
ATTN: R. POPE 
ATTN: T. MAZOLLA 
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SCIENCE APPLICATIONS INT'L CORPORATION 
P. O. BOX 1303 
MCLEAN, VA 22102 
ATTN: A. SARMA 
ATTN: C. AGRITELLIS 
ATTN: D. BACON 
ATTN: D. E. MAYS 
ATTN: J. COCKAYNE 
ATTN: J. MCGAGHAN 
ATTN: M. HALL 
ATTN: M. TURNER 
ATTN: N.AHMAD 
ATTN: P. L. R. MADALA 
ATTN: P. LEE 
ATTN: T. J. DUNN 

SCIENCE APPLICATIONS INT'L CORPORATION 
P. O. BOX 1303 
MCLEAN, VA 22102 
ATTN: T.WAIT 
ATTN: Y.JIN 
ATTN: Z. BOYBEYI 

THE AEROSPACE CORPORATION 

P. O. BOX 92957 

LOS ANGELES, CA 90009 

ATTN: DR M. PLONSKI, MS M8/719 

THE TITIAN CORPORATION 
TITAN RESEARCH & TECHNOLOGY DIVISION 
P. O. BOX 2229 
PRINCETON, NJ 08543-2229 
ATTN: I. SYKES 


TITAN CORPORATION (THE) 

TITAN RESEARCH & TECHNOLOGY DIVISION 
9410 TOPANGO CANYON BOULEVARD, 
SUITE 104 

CHATSWORTH, CA 91311-5771 
ATTN: R. ENGLAND 

TRW S. I. G. 

STRATEGIC SYSTEMS DIVISION 
P. O. BOX 1310 

SAN BERNARDINO, CA 92402-1310 
ATTN: N. LIPNER 


DEPARTMENT OF ENERGY 

UNIVERSITY OF CALIFORNIA 
LAWRENCE LIVERMORE NATIONAL 
LABORATORY 
P. O. BOX 808 

LIVERMORE, CA 94551-9900 
ATTN: L-030, A. KUHL 

LOS ALAMOS NATIONAL LABORATORY 

P. O. BOX 1663 

LOS ALAMOS, NM 87545 

ATTN: MS J514, A. S. MASON 
ATTN: MS670, J. NORMAN 
ATTN: R. W. WHITAKER, ESS-5MS/F665 
ATTN: WX-1, B.SHAFER 

DEPARTMENT OF THE AIR FORCE 

AF WEATHER TECHNICAL LIBRARY 
151 PATTONAVENUE 
ASHEVILLE, NC 28802 
ATTN: K. MARSHALL 


HEADQUARTERS 
AIR FORCE SPACE COMMAND 
DIRECTOR OF OPERATIONS 
150 VANDENBERG STREET 
PETERSON AFB, CO 80914-4120 
ATTN: AFSPC/DO 


AIR FORCE WEATHER AGENCY/DNXM 
106 PEACEKEEPER DRIVE, SUITE 2N3 
OFFUTT AFB, NE 68113-4039 
ATTN: MAJ R. LEFEVRE 

AIR UNIVERSITY LIBRARY 
600 CHENNAULT CIRCLE 
BLDG 1405-ROOM 160 
MAXWELL AFB, AL 36112-6424 
ATTN: AUL-LSE 


HQ USAF/XOWX 
1490 AIR FORCE PENTAGON 
WASHINGTON, DC 20330-1490 
ATTN: XOWX 


VISIDYNE, INC. 

P. O. BOX 1399 
GOLETA, CA 93116-1399 
ATTN: J. DEVORE 
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DEPARTMENT OF THE ARMY 


DEPARTMENT OF THE ARMY 
DEPUTY CHIEF OF STAFF FOR 
OPERATIONS AND PLANS 
PENTAGON 

WASHINGTON, DC 20310-0460 
ATTN: DAMO-NCZ 


US ARMY RESEARCH LABORATORY 
AMSRL-SL-CS E3331 
5101 HOADLEY ROAD 

ABERDEEN PROVING GROUND, MD 21010-5423 
ATTN: SLCBR-SS-T, TECHNICAL LIBRARY 


US ARMY RESEARCH LABORATORIES 
2800 POWDER MILL ROAD 
ADELPHI, MD 20783-1197 
ATTN: AMSRL-SL-CE 
ATTN: J. MARTIN 
ATTN: R. CIONCO 


COMMANDER 

WEST DESERT TEST CENTER 
ATTN: STEDP-WD-M 
DUGWAY, UT 84022-5000 
ATTN: C. BILTOFT 
ATTN: J. BOWERS 


DEPARTMENT OF THE NAVY 


DEPUTY CHIEF OF NAVAL OPERATIONS 
(PLANS, POLICY & OPERATIONS) 

2000 NAVY PENTAGON ROOM 
ROOM4E592, N3/N5 
WASHINGTON, DC 20350-2000 
ATTN: BRANCH HEAD, N514 

NAVAL RESEARCH LABORATORY 
4555 OVERLOOK AVENUE, SW 
WASHINGTON, DC 20375-5000 

ATTN: CODE 5227, RESEARCH REPORT 
ATTN: S. CHANG 

COMMANDER 

NAVAL SURFACE WARFARE CENTER 
DAHLGREN DIVISION 
17320 DAHLGREN ROAD 
DAHLGREN, VA 22448-5000 
ATTN: CODE BSI ,T. BAUER 




